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Abstract 

■  ■  - 

This  study  modifies,  the  thermodynamic  model  of  a 

r 

-  r” 

previously  existing  first-order  accurate  Flux-Dif ference<- 
Splitting  (FDS)  algorithm  for  planar,  supersonic  nozzles. 

The  thermodynamic  model  is  changed  from  a  calorically  and 
thermally  perfect  gas  to  a  thermally  perfect  (imperfect) 
gas,  where  the  flow  field  is  "frozen"  or  non  -  reacting .  The 
frozen  flow  and  imperfect  gas  assumptions  more  nearly 
approximate  the  real  behavior  of  a  fluid  in  supersonic 
propulsive  nozzles.  The  modified  code  can  now  account  for 
specific  heats  that  vary  as  a  function  of  temperature. 

Using  curve  fittings  of  JANAF  thermochemical  data,  the  code 
can  handle  -nine  gas  species,  as  well,  to  model  combustion 
products  entering  the  nozzle  inlet.  The  marching  scheme  is 
not  altered  in  order  to  retain  the  robustness  and  efficiency 
of  the  first-order  method. 

An  oblique  shock  reflection  study  is  done  to  validate 
the  improved  gas  model.  A  low  pressure,  low  temperature 
case  and  a  high  pressure,  high  temperature  case  are  run. 

For  the  first  case,  the  perfect  and  imperfect  models  are 
nearly  identical.  For  the  more  extreme  case,  pressure  for 
the  perfect  gas  is  9.4%  greater  than  the  exact  solution,  at 
the  upper  boundary,  across  the  shock.  An  interior  flow 

I  » 


X 


^nozzle  is  run  for  the  two'  cases,  with  air  as  the  working 
fluid.  Again,  the  two  models  give  identical  resi-'.ts  for  the 
low  pressure  case.  For  the  high  pressure  case,  integrated 
nozzle  thrust  for  the  imperfect  gas  model  is  16%  higher  than 
that  for  the  original  perfect  gas  model. 


IMPROVEMENT  OF  THE  THERMODYNAMIC  MODEL  FOR  A  FLUX- 


DIFFERENCE-SPLITTING  ALGORITHM  FOR  THE  COMPUTATION  OF 

HIGH-SPEED  FLOWS 

I .  Introduction 

1 . 1  Historical  Development 

From  the  earliest  days  of  space  exploration,  scientists 
and  engineers  have  explored  the  concept  of  a  single-stage- 
to-orbit  (SSTO)  flight  vehicle.  Such  a  design  could  takeoff 
conventionally  from  existing  runways,  enter  a  low-earth 
orbit,  reenter  the  earth's  atmosphere,  and,  finally,  land 
conventionally  on  a  runway.  To  do  so  v;ith  an  autonomous 
vehicle  stretches  the  current  limits  of  technology  in  the 
areas  of  structures,  materials,  aerodynamics,  and 
propulsion.  In  particular,  the  air-breathing,  hydrogen- 
fueled  propulsion  system  currently  proposed  must  perform 
over  the  wide-  range  of  flight  conditions  that  will  be 
encountered  during  a  typical  mission  profile.  This  places 
special  demands  on  researchers'  abilities  to  design  and  test 
propulsion  systems  before  actual  flight  occurs. 

The  tools  of  computational  fluid  dynamics  (CFD)  will 
provide  the  lion's  share  of  this  propulsion  testing 
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1. 

capability.  More  specifically,  CFD  will  provide  the  most 
economical  and,  at  times,  the  only  means  of  testing  nozzle 
designs  in  the  high  Mach  number  flight  regimes  (2:1). 

In  response  to  this  need  for *^mputational  design  tools 
to  model  engine  flow  fields,  a  computer  code  was  developed 
by  Doty  (3)  to  model  the  flow  of  a  calorically  and  thermally 
perfect  gas  as  it  applies  to  a  two-dimensional,  maximum 
thrust,  planar  nozzle.  Figure  1  shows  the  location  and  type 
of  geometry  for  a  superson'::,  planar  nozzle  on  a  NASP-type 
vehicle.  What  remains  to  be  determined,  though,  is  the 
effect  of  modeling  the  nozzle  flow  with  this  simplified 
perfect  gas  model.  At  the  high  temperatures  involved 
(approximately  T=  4000  K),  specific  heats  of  gas  are  not 
constant,  as  a  calorically  and  thermally  perfect  gas  implies 
(from  here  on  referred  to  as  a  perfect  gas),  but  vary 
significantly  as  a  function  of  temperature  (1:391).  This 
variation  can  have  a  measurable  effect  on  flow  field 
properties  and  thrust  calculations. 

1 . 2  Problem  Statement 

This  research  effort  incorporates  an  imperfect  gas 
model  (thermally  perfect)  into  a  previously  existing 
algorithm,  originally  based  on  a  perfect  gas,  to  more 
realistically  predict  the  performance  of  the  installed 
nozzle . 
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1 . 3  Objectives 


The  following  objectives  of  this  study  are  outlined 
below: 

i.  Become  familiar  with  the  numerical  schemes  and 
solution  procedures  used  in  the  reference  (3) 
nozzle  flow  program. 

ii.  The  FDS  method  requires  the  solution  of 

discontinuities  in  the  flow  field,  also  known  as 
solving  the  Riemann  problem.  With  the  aid  of 
curve-fitted  chemical  composition  data  from  an 
existing  software  package,  implement  varying 
specific  heats  of  gas  for  hydrogen/air  combustion 
products  into  the  solution  of  the  Riemann  problem. 
The  basic  numerical  solution  procedure  will  remain 
unchanged. 

iii.  Validate  the  modified  program  against  an  exact, 
imperfect,  oblique  shock  reflection  test  case. 

iv.  Compare  thrust  calculations  of  the  perfect  gas 
case  to  the  imperfect  gas  case.  A  determination 
will  then  be  made  as  to  whether  the  improved 
model  provides  enough  increase  in  accuracy  to 
justify  the  added  computational  time. 

1 . 4  Assumptions 

The  primary  assumption  of  this  study  is  that  the  gas 
obeys  the  thermally  perfect  equation  of  state: 
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P  =  pRT 


This  relation  holds  for  either  a  perfect  or  imperfect  gas. 
Thermochemically ,  the  gas  flow  composition  will  be  fixed  at 
some  initia.l  value  before  solution  of  the  nozzle  flow  field 
begins.  This  would  be  the  case  if  chemical  reactions  could 
be  inhibited  while  still  allowing  collisions  to  occur.  E’low 
properties  such  as  pressure,  temperature,  density,  and 
specific  heats  may  still  vary.  This  inhibited-reaction  flow 
is  referred  to  as  "frozen"  flow  (9:191-192).  This 
assumption  is  justified  due  to  the  extremely  low  pressures 
and  high  speeds  that  gas  molecules  encounter  at  high 
altitudes.  Low  pressure  at  higher  altitudes  contribute  to 
large  molecular  mean  free  paths  (distance  between 
collisions),  lessening  the  chances  of  collisions  resulting 
in  reactions.  Additionally,  molecules  traveling  at  high 
speeds  relative  to  some  nozzle  length  scale  will  have 
extremely  small  nozzle  residence  times  and  will,  therefore, 
have  little  time  to  chemically  react  with  each  other. 

1 . 5  Scope 

This  study  will  deal  strictly  with  modifying  the  code's 
current  thermodynamic  model.  The  numerical  scheme  employed 
will  remain  intact.  While  the  program  has  the  capability  to 
design  maximum  thrust  nozzles  and  predict  nozzle  performance 
over  a  wide  range  of  nozzle  parameters,  a  standard  nozzle 
configuration  and  specified  flight  condition  will  be  used 
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during  thermodynamic  model  comparisons . 


1 . 6  Background 

The  recency  of  work  done  in  reference  (3)  meant  that 
the  chances  of  finding  a  sudden  rise  in  literature  on  the 
FDS  method  and  it's  application  of  the  Riemann  problem  were 
not  high.  A  thorough  search  of  the  work  done  regarding  CFD 
algorithms  and  their  attendant  thermodynamic  models  was 
performed.  No  new  work  has  been  done  in  the  area  of 
internal  flow  (nozzles)  using  a  first-order  FDS  method  with 
local  solution  of  the  Riemann  problem  to  handle  flow 
discontinuities.  By  default,  the  improvement  of  its 
thermodynamic  model  would  be  an  original  effort.  But  many 
approaches  to  thermodynamic  modeling  were  surveyed  whether 
or  not  they  involved  FDS. 

Computational  methods  usually  use  the  perfect  gas  model 
for  validation  of  the  scheme.  Subsequent  efforts  to  improve 
the  gas  model  usually  go  beyond  an  imperfect  gas  with  frozen 
flow  and  straight  to  one  that  can  handle  finite-rate 
chemistry  and  vibratjonal  nonequilibrium  (8:1).  Many  can 
account  for  dissociating  gases.  While  improved  correlation 
to  physics  is  the  goal,  there  is  an  accompanying  increase  in 
CPU  time  which  can  be  very  expensive  if  the  code  must  be  run 
on  supercomputers. 
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Nozzle 


Figure  1  Hypersonic  vehicle  and  planar  nozzle  (3:5) 
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I I .  Governing  Equations  and  Coordinate  Transformaticn 


2 . 1  Governing  Equations 

The  governing  equations  of  motion  for  a  planar,  steady, 
adiabatic,  inviscid  flow  of  a  compressible  fluid  with  no 
external  work  or  body  forces  are  the  Euler  equations, 
written  here  in  divergence  vector  form  as : 

=  o  (11 

bx  ay 

where  the  E  and  F  vectors  are  given  in  terms  of  the 
conservation  variables  as: 


pu 

_  pu^+P 
puv 

u(pe+P) , 


pv 
pvu 
pv^+P 
v(pe+P) , 


The  first  term  of  each  vector  represents  continuity,  the 
second  and  third  are  the  x  and  y  components  of  momentum, 
respectively,  and  the  fourth  is  the  energy  f^quation.  The 
governing  vector  equation  holds  true  for  both  a  perfect  and 
imperfect  gas. 


2 . 2  Thermodynamic  Model 


The  original  code  uses  a  thermodynamic  model  based  on 
a  thermally  and  calorically  perfect  gas.  As  such,  the 
thermal  equation  of  state  for  this  perfect  gas  assumption 
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P  =  pRT 


(3) 


The  total  internal  energy  for  a  perfect  gas  is: 

pe  =  +  -|p(u2  +  v^)  (4) 

The  thermally  perfect  and  calorically  imperfect  gas 
(imperfect  gas)  assumptions  and  equations  are  presented  in 
Appendix  A.  The  thermal  equation  of  state,  eq.  (3),  remains 
valid  for  an  imperfect  gas.  But  the  total  specific  internal 
energy  shown  uses  perfect  gas  assumptions  to  obtain  the 
first  term  on  the  right  hand  side.  The  total  internal 
energy  for  an  imperfect  gas  is: 

pe  *  ph  -  P  +  —p{u^  +  v^)  (5) 

Static  enthalpy,  h,  is  obtained  from  a  least  squares  curve 
fit  of  thermochemical  data  from  JANAF  tables  (4:32).  The 
enthalpy  of  the  mixture  is  given  as  (10:50): 

h  =  hg  +  f^CpdC  (6) 

where  the  specific  heat,  Cp,  given  as  a  function  of 
temperature  and  the  least  squares  coefficients,  is  of  the 
form  (4:32)  : 

Cp  =  (a  +  bT  *  cT^  *  dT^  +  eT^)Rg^^  (7) 
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2 . 3  Coordinate  Transformation 

The  numerical  solution  ic  accomplished  in  the 
computational  plane  and  requires  that  the  governing 
equations,  eq .  (1),  be  transformed  to  this  space.  The 
physical  coordinates  x  and  y  are  mapped  to  the  computational 
coordinates  as  follows  (3:142): 

C  =  X  T)  =  T)  (x,y)  (8) 

The  governing  equations,  eq.  (1),  are  transformed  in  uniform 
computational  space  as  (3:8): 


d{E)  _  d{E)  d{F) 

aC  3ti  dr] 


(9) 


For  further  details  on  the  coordinate  transformation,  see 
(3:142) . 


Ill .  Numerical  Algorithm 


3 . 1  Introduction 

The  Flux-Difference-Splitting  (FDS)  method  is  a 
technique  that  requires  the  solution  of  the  Riemann  problem 
to  locally  iTiodel  the  correct  physics  of  the  flov/  and 
globally  handle  complex  flow  interactions.  Figure  2 
illustrates  the  Riemann  problem.  The  Riemann  problem  is  the 
solution  of  flo..'  discontinuities,  and  the  waves  that  result 
propagate  information  in  specific  directions  based  on  the 
wave  angle  (3:9).  Fluxes  across  waves  can  be  calculated 
using  known  values  in  region  (6)  and  (0)  and  the  solution 
aft  of  the  discontinuity.  These  fluxes  are  then  "split" 
based  on  the  propagation  direction  (flow  angle),  in  effect 
weighting  the  flow  information  in  the  physically  correct 
direction  to  the  next  solution  plane,  i+1.  Figure  3 
illustrates  the  computational  space. 

A  first-order  accurate  marching  scheme  is  used  with  the 
FDS  algorithm  for  its  monotonic,  or  nonoscillatory ,  behavior 
and  robustness.  As  described  by  van  Leer  and  reported  by 
Doty,  in  regions  of  strong  gradients,  this  is  very  important 
(3:10).  For  further  details  on  the  first-order  accurate  FDS 
scheme,  refer  to  Appendix  C  and  reference  (3).  This 
combination  of  proper  local  physics  and  first-order  behavior 
allows  strong  gradients  (shock  waves  and  contact  surfaces) 
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to  be  handled  very  accurately.*  (3:2). 

3 . 2  The  Riemann  Problem 

Figure  2  illustrates  the  setup  of  the  Riemann  problem. 
As  described  by  Doty,  Godunov  proposed  that  the  general  flow 
property,  W,  be  modeled  as  a  series  of  uniform  flow  regions 
(3:11).  The  discontinuity  is  assumed  to  be  located  halfway 
between  node  points  j  and  j+1  at  Riemann  node  j+1/2.  Waves 
(1)  and  (3)  can  be  any  combination  of  expansion  or 
compression  waves  (or  neither,  in  which  case  there  is  no 
turning  of  the  flow) .  Wave  (2)  is  a  contact  surface  that 
separates  regions  (2)  and  (4).  Properties  are  known  at  the 
nodes  j  and  j+1,  which  correspond  to  regions  (0)  and  (6), 
respectively  (3:11). 

3 . 3  Solution  of  the  Riemann  'Problem 

The  problem  consists  of  solving  for  the  pressure, 
density,  and  velocity  components  in  regions  (2)  and  (4). 

This  is  done  by  one  of  three  methods.  The  first  is  an  exact 
solution  that  requires  an  iterative  procedure  to  solve  non¬ 
linear,  coupled,  non- isentropic  relations  across  compression 
waves  (shock  waves)  and  Prandtl-Meyer  (simple)  expansion 
waves.  Besides  the  iteration  of  the  non-linear  equations, 
the  solution  of  pressure  in  regions  (2)  and  (4)  must  be 
solved  iteratively  because  contact  surfaces  can  not  support 
normal  pressure  jumps  and  the  pressures  must,  therefore. 
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match  across  the  contact  surface.  The  second  method  is  the 
exact-approximate  solution.  This  is  very  similar  to  the 
exact  solution  except  compression  waves  are  treated  as 
isentropic  compressions.  Again,  non-linear  equations,  this 
time  due  to  Prandtl-Meyer  relations,  must  be  solved 
iteratively  and  contact  surface  requirements  must  be  met. 

For  details  on  these  exact  methods,  see  reference  (3).  The 
last  is  the  linearized-approximate  method  and  is  the  focus 
of  this  study. 

3.3.1  Linearized-Approximate  Method 
This  approach  eliminates  the  iteration  needed  in  the 
previous  methods  by  linearizing  the  Prandtl-Meyer  relations. 
In  addition,  compression  and  expansion  waves  are  assumed  to 
be  isentropic.  See  Appendix  B  for  details  on  this  method. 

It  starts  with  the  differential  form  of  the  compatibility 
relations : 

dP  Jb  pV^dfd  =  0 

where  the  velocity  magnitude  and  flow  angle  are  defined  as 
follows : 

(11) 

6  =  tan'*  ( v/u)  ( ) 

Through  manipulation,  the  linearized  equation,  for  the  case 
where  wave  (3)  is  a  compression  wave  and  wave  (1)  is  an 
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expansion  wave,  becomes,  respectively: 

[ln(P)]4  +  =  [ln(P)]6  + 

[ln(P)]2  -  (ZqJoz  = 

where  the  linearized  coefficient  is: 

z  =  lY^Va^)  ,15, 

Recall  that  the  pressures  and  flow  slope  must  match  across 
the  contact  surface: 

P4  =  P2  (16) 

Of  =  O2  (17  ) 

Equations  (13),  (14),  (16),  and  (17)  represent  four 
equations  and  four  unknowns  and  may  be  solved  simultaneously 
to  yield  values  for  and  o,. 

3, 3. 1.1  Perfect  Gas. 

Solving  the  simultaneous  equations  for  pressure  in 
region  (4),  P4,  yields: 

P4  =  exp(  (ln(P)  ]g  +  Zg  (o^-Og) )  (18) 

All  terms  on  the  right  hand  side  are  known.  Properties  in 
region  (6)  are  known  since  these  are  initial  values  at  plane 
i.  The  Zg  coefficient  can  be  "lagged"  in  the  known  region 
(6)  or  updated  in  region  (4)  to  provide  a  more  accurate 
value  of  the  coefficient.  Recall  that  the  slope  of  the  flow 
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in  region  (4),  O4,  is  found  from  the  simultaneous  solution 
of  eqs.  (13),  (14),  (16),  and  (17)  and  is  given  as: 


[ln(P)  ]  Q  -  [ln(P)  ]g+(Ze)  Og  +  (Zq)  Op 


(19) 


The  details  of  this  derivation  appear  in  Appendix  B  and 
reference  (3). 

With  pressure  known  aft  of  the  wave,  isentropic 
relations  such  as: 


Pi  = 
Pe 


(20) 


34  =  (21) 

can  be  used  to  obtain  density  and  speed  of  sound  in  region 
(4).  For  adiabatic  flow,  conservation  of  stagnation 
enthalpy  (hf^  =  across  wave  3  is  used  to  determine  Mach 

number  and  the  velocity  components  in  regions  (4).  But 
these  rely  on  a  constant  specific  heat  ratio,  y,  so  another 
approach  must  be  taken  for  an  imperfect  gas.  See  reference 

(3)  for  details  of  the  perfect  gas  Riemann  problem  solution. 

3. 3. 1.2  Imperfect  Gas. 

All  the  above  manipulations  and  equations  from  eq.  (10) 
to  eq.  (18)  are  valid  for  an  imperfect  gas  as  well. 
Therefore,  to  obtain  the  other  property  variables  in  region 

(4) ,  a  relationship  is  needed  that  can  utilize  pressure  as 
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the  only  known  value  in  region  (4).  The  following  equations 
show  how  this  requirement  can  be  met  when  the  waves  are 
approximated  as  isentropic  waves.  Across  the  isentropic 
wave  3,  going  from  region  (6)  to  region  (4),  the  change  in 
entropy  across  the  shock  is  zero  (10:58); 

S^-Sg  =0  (22) 

Integrating  the  differential  entropy  equation  yields: 


Combining  theses  two  relations  results  in  the  entropy 
parameter  in  region  (4), 

Kin  (24) 

The  pressure  in  region  (6),  P^,  is  known  and  P,  was 
calculated  from  eq.  (18).  The  entropy  parameter  in  region 
(6),  (|>6,  is  calculated  using  the  least  squares  coefficients 
and  the  known  temperature  in  region  (6)  (4:32): 

+  alnTg  >  bT^  *  ^  |r*"]  R  (25) 

With  a  value  for  <|>4  given  by  eq.  {2A) ,  the  following 
equation  for  «(>4f  a  function  of  temperature  in  region  (4), 
is  given  in  least  squares  coefficient  form  (4:32): 
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(26) 


<f)4  =  [(j)„  +  alnT^  +  bT^  +  ^  ^ 


This  equation  can  be  iteratively  solved  for  temperature  in 
region  (4),  T^,  and  the  equation  of  state  then  determines 
density . 

For  adiabatic  flow,  the  imperfect  case  also  relies  on 
conservation  of  stagnation  enthalpy  (h^g  =  ht4)  to  determine 
velocity  components  and  Mach  number,  but  static  enthalpy 
must  be  determined  as  a  function  of  temperature  in  order  to 
find  the  x  component  of  velocity,  u.  Total  stagnation 
enthalpy  in  region  (6)  and  (4)  is: 


(27) 


The  local  slope  of  the  flow  in  region  (4),  Oj,  is  defined 


O4  =  V4/U4  (28) 

Substituting  this  equation  into  eq .  (27)  and  solving  for  the 
X  velocity  component,  u^,  yields: 


I+O4 


where  static  enthalpy  in  region  (4),  h^,  can  be  found  as  a 
function  of  T^  and  the  least  squares  coefficients.  More 
complete  details  of  tiie  imperfect  gas  Riemann  problem 
solution  are  contained  in  Appendix  B. 


3 . 4  First-Order  Accurate  Interior  Point 


The  initial  values  are  known  at  plane  i  and  the 
solution  is  desired  at  the  next  plane,  i+1.  Figure  3  shows 
the  stencil  in  the  solution  plane.  The  first-order  accurate 
FDS  method  uses  positively  biased  information  from  Riemann 
node  j-1/2  and  negatively  biased  information  from  Riemann 
node  j+1/2  to  integrate  the  solution  from  node  (i,j)  to 
(i,j+l).  To  march  the  solution  in  the  x  direction,  a 
combination  of  finite  difference  and  flux  difference 
approximations  are  used  for  the  partial  derivatives  of  the  E 
and  F  vectors  in  the  governing  equation.  This  is 
accomplished  by  computing  the  Riemann  fluxes  and  then 
differencing  these  fluxes.  These  flux  differences  are  then 
propagated  or  "split"  in  a  particular  direction  based  on  the 
local  slope  of  the  flow.  It  should  be  noted  that,  after 
solving  the  Riemann  problem  by  the  new  method  described  in 
the  previous  section,  this  FDS  method  remains  unchanged. 
Details  of  the  solution  procedure  are  in  Appendix  C. 

The  first-order  accurate  FDS  approximation  to  the 
transformed  governing  equation  is  (3:15): 

[dE;.u2*dEj,,/^]-M:r\y  > 

The  first-order  accurate  FDS  boundary  point  calculations  are 
presented  in  reference  (3). 
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3 . 5  Decoding  the  Solution 

With  the  solution  vector  integrated  from  plane  i  to 
i+1,  the  E  and  F  flux  vectors,  in  conservative  form,  are 
decoded  or  reduced  to  the  primitive  variable  form  necessary 
for  the  solution  of  the  Riemann  problem  at  the  next  plane. 
Details  are  presented  in  Appendix  D.  The  E  vector,  in 
conservative  form,  was  presented  in  Chapter  2  and  is 
repeated  here  for  convenience: 

pu 

E=  (31) 

puv 

u(pe+P). 

This  can  be  written  in  terms  of  the  E  vector  components, 
which  are  known  from  the  marching  scheme: 


El  =  pu 

(32) 

£2  =  pu^  +  P 

(33) 

E3  =  puv 

(34) 

E4  =  u(pe  +  P) 

(35) 

As  was  mentioned  in  Chapter  2,  the  E4  vector  differs  in 
total  internal  energy,  pe,  for  the  two  thermodynamic  models 
For  a  thermally  perfect  gas,  total  internal  energy  is  given 
as : 
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(36) 


1 

pe  =  pu  +  — p  (u^  +  v^) 

where,  for  a  perfect  gas,  specific  internal  energy,  6,  is: 

u  =  c^T  (37) 

and  for  an  imperfect  gas,  specific  internal  energy  is: 

n  =  h  -  RT  (38) 

After  substituting  eq.  (37)  into  eq.  (36)  and 
performing  several  manipulations,  decoding  of  the  perfect 
gas  solution  provides  a  quadratic  equation  for  the  x 
component  of  velocity,  u,  which  is  a  function  of  the  E 
vector  components  and  y  (3:222): 

Y_tl  (El)  -  [-1-  {E2)  1  m-  [  (g4)  -1  -7^1  '  ° 

2  (y-l)  [  Y-1  2  {ED 

Since  this  equation  relies  on  a  constant  y,  another 
approach  is  required  for  an  imperfect  gas.  Details  of  the 
imperfect  gas  decoding  procedure  are  presented  in  Appendix 
D.  Begin  by  substituting  eq.  (38)  into  eq .  (36).  Substitute 
the  resulting  eq.  (36)  into  the  E4  vector,  eq .  (35).  After 
manipulating  .and  substituting  eqs.  (32)  and  (34)  into  eq. 
(35),  the  E4  vector  component  becomes: 

E4  =  hEl  +  -Elu^  + 

2  2  El 

This  equation  has  two  unknowns:  static  enthalpy,  h,  and  the 
velocity  component,  u.  Static  enthalpy  is  a  function  of 
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temperature  and  the  least  squares  coefficients.  By 
substituting  the  thermal  equation  of  state  and  the  E2 
component  into  the  El  component,  the  x  component  of 
velocity,  u,  can  also  be  cast  as  a  quadratic  in  terms  of  the 
unknown  temperature : 

(El)  (E2)  U  (El)  RT  =  0  (41) 

An  initial  guess  for  temperature  can  now  be  put  into  eq. 

(41)  and  static  enthalpy,  h.  If  the  E4  component  is  not 
satisfied,  the  process  is  repeated.  An  iterative  process, 
such  as  the  secant  method,  can  be  used  to  find  temperature. 
The  calculation  of  the  other  primitive  variables  is  as 
follows : 


E2  ±  \/Ts2r^4TsirrrEiyRiT 

(42) 

2{E1) 

(43) 

El 

P  =  E2  -  (El)  U 

(44) 

p  =  Il 
u 

(45) 

In  summary,  the  Riemann  problem  is  solved  at  all  the 
half  node  locations  at  plane  i  by  one  of  the  three  solution 
methods  described.  The  fluxes  are  differenced  and  split  and 
then  marched  by  the  first-order  accurate  method  shown.  At 
plane  i+1,  the  solution  vector,  E,  is  decoded  to  start  the 
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process  again. 

This  approach  differs  from  the  more  familiar  marching 
algorithms.  Where  decoding  the  E  vector  into  primitive 
variables  is  optional  for  methods  such  as  the  Roe  scheme, 
the  Riemann  problem  requires  it.  Again,  the  purpose  of  this 
is  to  model,  more  directly,  the  correct  physics  of  the  local 
flow. 
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(d^d_F)]+1/2 

(dE.dF)]_i/2 


IV.  Validation  Studies 


4 . 1  Introduction 

The  imperfect  gas  model  was  validated  by  comparing 
results  to  exact  solutions.  There  is  no  exact  solution  for 
the  flow  field  in  a  hypersonic  nozzle.  But  there  are  exact 
solutions  available  for  various  flow  geometries.  One  is  the 
oblique  shock  reflection.  Reference  (10)  has  a  parallel 
development  structure  that  discusses  perfect  gas  relations 
and  the  corresponding  imperfect  gas  case  for  normal  shocks, 
oblique  shocks,  and  Prandtl-Meyer  waves.  The  oblique  shock 
case  v/as  chosen  since  the  large  jump  in  properties  across 
such  a  shock  are  a  good  test  of  the  numerical  method.  For 
details  on  validation  of  the  perfect  gas  algorithm,  using 
all  three  Riemann  solvers,  for  supersonic  source  flov;  and 
Prandtl-Meyer  flow,  see  reference  (3). 

4 . 2  Oblique  Shock  Reflection  Study 

The  geometry  for  the  oblique  shock  reflection  study  is 
shown  in  Figure  4.  The  height  of  the  channel  is  0.0254 
meters  before  the  ramp  and  the  length  of  the  channel  is  0.11 
meters.  Computations  begin  at  point  0  and  proceed  to  a 
location  past  the  second  dashed  line.  Uniform  flow  travels 
left  to  right  and  encounters  an  abrupt  turn  at  the  10°  ramp, 
depicted  as  The  resulting  incident  shock  impinges  the 
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ceiling  of  the  channel  and  a  reflected  shock  is  produced. 

The  region  used  for  validation  is  between  the  dashed  lines. 

A  relatively  fine  grid  of  51  node  points,  in  the  y 
direction,  is  used  for  all  cases.  Grid  spacing  in  the  x 
direction  is  based  on  the  Courant-Friedrichs -Lewy  (CFL) 
stability  criteria  (3:144).  Details  on  grid  generation  can 
be  found  in  reference  (3). 

Two  initial  conditions  were  run  for  comparison  of  the 
two  thermodynamic  models.  Input  for  the  initial  value  (IV) 
line  consists  of  the  x  and  y  location  of  the  IV  line,  static 
pressure,  static  temperature,  Mach  number  and  flow  angle. 

The  data  for  initial  condition  #2,  IC  2,  was  obtained  from  a 
cycle  analysis  code  that  outputs  species  concentrations  and 
flow  properties  (6).  Table  I  shows  the  tv/o  initial 
conditions . 

For  the  validation  studies,  air  was  the  working  fluid 
in  both  thermodynamic  models.  The  perfect  gas  model  is  very 
limited  in  its  ability  to  model  different  fluids.  This  is 
accomplished  in  the  input  deck  by  specifying  the  specific 
heat  ratio,  7,  and  the  gas  constant,  where  R,j«k= 

R,„uv/mole-wt .  The  perfect  gas  model  sets  7=  1.4  and  R, 

287.0,  for  the  validation  study.  The  imperfect  gas  model 
uses  a  more  chemically  correct  air  composition  in  the 
familiar  mass  fractional  quantities: 


2 


^2 

Oz 

Ar 


.7556 

.2316 

.0128 


mix 

^9o 


^9mix 

^9Ar 

^9„ix 


The  gas  constant  is  a  function  of  the  molecular  weight  of 
the  gas  and  is,  therefore,  the  same  value  as  that  for  the 
perfect  gas,  but  the  specific  heat  ratio  is  a  function  of 
temperature.  Recall  that  the  fluid  composition  is  unchanged 
or  "frozen"  throughout  the  flow  field. 

Two  areas  in  the  flov;  were  studied.  The  first  was  the 
ceiling  of  the  channel  or  the  top  most  node.  Here,  it  could 
be  seen  whether  the  full  pressure  rise  on  both  sides  of  the 
shock  impingement  point  would  be  captured  and  located 
correctly.  The  second  was  at  interior ’ node  j=  40  (node 
counting  ascends  from  the  floor  to  the  ceiling)  .  This 
interior  node  location  will  shov/  the  capture  and  location  of 
the  incident  and  reflected  shock  waves. 

The  exact  solution  of  an  oblique  shock  wave  used  is 
based  on  an  imperfect  gas.  This  solution  requires  that  a 
double  iteration  be  performed  on  values  of  the  shock  wave 
angle,  e,  and  density,  p.  Figure  5  illustrates  the  oblique 
shock  wave  geometry.  Details  of  this  method  are  found  in 
reference  10  (10:369).  A  short  code  v/as  developed  to 
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quickly  find  properties  aft  of  an  imperfect  oblique  shock  to 
use  as  input  for  the  initial  value  line  of  the  code.  This 
program  is  located  in  Appendix  E. 

Figure  6  shows  the  results  of  IC  1  for  the  ceiling 
boundary  condition.  The  two  models  show  no  discernable 
difference  in  pressure  between  the  two  solutions.  Both 
models  correctly  capture  the  magnitude  of  the  pressure  jump 
and  locate  the  shock  well.  The  agreement  is  expected  given 
the  low  pressure  and  low  temperature  initial  conditions.  At 
temperatures  less  than  T=  600K,  y  remains  relatively 
constant  and  the  two  models  should  behave  similarly  (1:441). 
This  holds  true  for  T=  578  K,  where,  in  region  3,  y=  1.4  for 
the  perfect  gas  and  y=  1.39  for  the  imperfect  gas.  Looking 
at  Figure  7  for  the  interior  node,  again,  the  pressure 
magnitudes  and  locations  of  the  shocks  are  computed 
correctly . 

Figure  8  shows  the  results  for  IC  2  at  the  ceiling. 

For  this  boundary,  the  imperfect  gas  correctly  captures  the 
magnitude  and  location  of  the  pressure  jump,  but  the  perfect 
gas  misses  both  the  location  and  the  pressure  jump  across 
the  shock  by  a  significant  amount:  a  9.4%  pressure 
difference.  The  difference  between  the  thermodynamic  models 
can  be  attributed  to  the  greater  effect  that  the  higher 
temperature  has  on  y.  For  a  temperature,  T=  4278  K,  in 
region  3,  y=  1.276  for  the  imperfect  gas  and  y=  1.4  for  the 
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perfect  gas.  This  is  due  to  the  perfect  gas  having  a  higher 
percentage  of  energy  conversion,  across  the  shock,  show  up 
as  translational  molecular  energy  than  would  occur  for  an 
imperfect  gas.  Since  temperature  is  a  direct  measurement  of 
the  translational  energy  mode,  the  perfect  gas  model  has  a 
greater  magnitude  change  in  pressure  and  temperature  across 
the  shock.  For  details  on  the  mechanisms  of  high- 
temperature  effects  for  an  imperfect  gas,  see  Appendix  A. 

The  perfect  gas  model  misses  the  location  of  the  shock 
because  it  computes  a  shock  wave  angle,  £,  that  is  greater 
than  that  for  the  imperfect  gas.  The  shock  wave  angle  is, 
therefore,  located  a  shorter  distance  in  the  x  direction 
along  the  upper  wall. 

Figure  9  shows  the  results  for  the  interior  node  for  IC 
2.  Again,  the  imperfect  model  correctly  captures  the  shock 
location  and  pressure  magnitude  rise  for  both  waves.  The 
pressure  difference  between  the  exact  solution  and  the 
perfect  gas,  in  region  2,  is  5.3%  and  in  region  3,  9.4%. 

These  results  also  show  the  superb  behavior  of  the 
first-order  FDS  method.  The  most  obvious  characteristic  ij 
the  monotonic  behavior  of  the  method.  No  overshoot  of  th« 
pressure  magnitude  occurs  as  does  with  some  second-order 
accurate  methods  that  have  not  been  modified  to  detect 
gradients  and  ’•educe  oscillations  (3:24). 

The  preceeding  validation  studies  have  shown  that  the 
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improved  thermodynamic  model  has  been  correctly  implemented 
into  the  FDS  method  and  provides  consistent  results.  More 
physically  correct  solutions  are  now  possible  for  the 
conditions  likely  to  be  encountered  in  SCRAMJET  nozzles. 
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Figure  6  Pressure  along  upper  wall  for  IC  1 


Figure  7  Pressure  at  interior  node  j=40  for  IC  1 
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X- local  ion  (m) 

Figure  8  Pressure  along  upper  wall  for  IC  2 


X- location  (m) 

Figure  9  Pressure  at  interior  node  j=  40  for  IC  2 
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Table  1.  Initial  conditions  for  oblique  shock  reflection 
study 


IC  1 

IC  2 

ramp  angle  (e„.ui) 

0 

O 

O 

0 

Mach  Number 

6.0 

3.465 

pressure  (N/m^) 

1696 . 4 

200703 

temperature  (K) 

273.23 

3079.1 
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V .  Results  and  Discussion 


5 . 1  Introduction 

Analysis  of  different  flow  conditions  for  a  supersonic 
planar  nozzle  is  presented  in  this  chapter.  The  first  case 
is  for  the  internal  flow  of  the  supersonic  nozzle.  The 
second  case  involves  an  interior/exterior  flow  to  simulate 
actual  flight  conditions.  Static  pressure,  static 
temperature,  and  thrust  are  plotted  as  a  function  of  the  x 
location  along  the’  nozzle  wall. 

Doty's  code  (reference  3)  can  run  a  variety  of  nozzle 
geometries  ranging  from  a  straight  wall  nozzle  to  a  skewed- 
parabola  nozzle  wall.  The  lower  cowl  can  also  be  angled  at 
various  "hinge"  points  to  further  aid  nozzle  optimization. 
For  both  flow  configurations,  a  straight  cowl  and  a 
parabolic  nozzle  wall  are  used. 

5 . 2  Internal  Nozzle  Analysis 

Figure  10  illustrates  the  internal  nozzle  geometry. 

For  the  internal  flow  nozzle,  the  cowl  is  extended  the 
length  of  the  nozzle,  in  effect,  splitting  the  flow.  An 
attachment  angle,  ©p,  must  be  specified  for  all  nozzles; 

0B=  25°  is  arbitrarily  chosen  for  this  analysis.  Two  sets 
of  initial  conditions  were  run.  The  first  set  represents  a 
less  extreme  temperature  and  pressure.  The  second 
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represents  conditions  consistent  with  hypersonic  flight  at  a 
freestream  Mach  number  of  MU=  15.  The  later  conditions  were 
obtained  from  a  cycle  analysis  code  (6).  The  modified 
imperfect  gas  code  was  designed  to  simulate  nine  gas  species 
for  a  given  flow  field  but  the  three  constituents  of  air 
(nitrogen,  oxygen,  and  argon),  in  the  correct  mass  ratio, 
were  used  for  comparison  purposes.  Nozzle  inlet  conditions 
are  listed  in  Table  2. 

The  pressure,  temperature,  and  thrust  plots  for  IC  1 
are  shown  in  Figures  11,  12,  and  13,  respectively.  For  this 
low-pressure,  low- temperature  case,  the  two  gas  models  give 
nearly  identical  results.  This  behavior  is  expected  since 
the  variation  of  specific  heats  with  temperature  is 
negligible  at  these  conditions.  This  is  consistent  with  the 
IC  1  validation  study  case.  IC  1,  therefore,  proves  that 
the  imperfect  gas  model  is  consistent  with  the  proven 
perfect  gas  model  from  reference  (3).  The  plots  of  pressure 
and  temperature  illustrate  particular  characteristics  of  the 
nozzle.  The  large  drop-off  in  temperature  and  pressure,  in 
Figures  11  and  12,  respectively,  is  a  result  of  the  rapidly 
expanding  flow  around  the  nozzle  attachment  radius  (AB  in 
Figure  10)  .  The  subsequent  small  rise  in  pressure  and 
temperature,  is  a  result  of  recompression  of  the  flow  in  the 
parabolic  region  of  the  nozzle  (BC  in  Figure  10).  Figure  13 
shows  the  integrated  thrust  along  the  nozzle  wall,  where 
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both  models  predict  the  same  result. 

Pressure,  temperature,  and  thrust  are  plotted  for  IC  2 
in  Figures  14,  15,  and  16,  respectively.  For  this  case, 
noticeable  diff.erences  are  seen  between  the  thermodynamic 
models . 

Temperature  shows  a  greater  difference  in  Figure  15. 
Specific  heat  ratio,  7,  for  the  perfect  gas  is  held  constant 
at  7=  1.4  while  7  varies  with  temperature  for  the  imperfect 
gas.  Notice  also  that  temperature  drops  a  larger  amount  for 
the  perfect  gas.  Referring  to  Appendix  A,  energy 
conversion,  across  expansion  or  compression  waves,  is 
distributed  over  different  molecular  energy  modes,  with  a 
perfect  gas  receiving  a  greater  percentage  of  translational 
energy  (ie.,  temperature)  than  an  imperfect  gas  does.  That 
is  why  temperature  and  pressure  have  a  greater  decrease  in 
the  expansion  region  of  the  nozzle  for  the  perfect  gas 
compared  to  the  imperfect  gas.  Pressure  shows  the  least 
difference  in  Figure  14.  Recall  that  the  calculation  of 
pressure  for  the  linear-approximate  method  was  nearly  the 
same  for  both  gas  models.  The  only  difference  was  in  the 
calculation  of  the  z  coefficient  in  eq .  (15),  where  7  is  a 
function  of  temperature.  The  coupling  of  temperature  and 
pressure  in  the  thermal  equation  of  state  results  in  a 
similar  trend  for  pressure,  but  to  a  lesser  degree. 
Consequently,  pressure  differences  will  not  manifest 
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themselves  to  a  large  extent. 

A  smaller  magnitude  drop  in  pressure,  for  the  imperfect 
gas,  results  the  in  greater  thrust  shown  in  Figure  16. 

Recall,  in  the  shock  validation  study,  that  the  perfect  gas 
model  overpredicted  the  pressure  magnitude  rise  for  IC  2. 
This  trend  continues  for  the  interior  nozzle,  except 
pressure  is  falling  instead  of  rising.  The  main  point  is 
the  overprediction  of  magnitude  changes  in  properties.  In 
this  case,  integrated  thrust  along  the  nozzle  wall  for  the 
imperfect  gas  is  sixteen  percent  greater  at  the  nozzle  exit. 
This  is  a  significant  difference  and  illustrates  the  effect 
of  accounting  for  property  variations  as  a  function  of 
temperature. 

5 . 3  Internal/External  Nozzle  Analysis 

The  internal/external  flow  configuration  was  run  v/ith 
air  for  both  the  interior  and  exterior  regions.  Although 
the  multiple  species  from  combustion  products  of  H^air  can 
be  simulated,  the  variable  specific  heat  trends  for  air  are 
similar.  Figure  17  illustrates  the  nozzle  flow  and 
geometry.  Freestream  conditions  consisted  of  an  altitude  of 
129,900  m,  pressure  of  303.36  N/m^,  temperature  of  250  K, 
and  Mach  number  of  15.0.  Initial  conditions  for  the 
exterior  flow  were  obtained  by  sending  air  across  a  10° 
ramp,  simulating  flight  at  some  angle  of  attack.  This  is  a 
simplification  of  the  forebody  compression  that  would 
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actually  occur.  These  initial  conditions  will  make  up  the 
properties  in  region  2  of  Figure  17.  Perfect  ani  imperfect 
gas  relations  were  used  for  computing  external  initial 
conditions  for  the  respective  thermodynamic  models .  Table  3 
shows  the  exterior  initial  conditions.  Table  4  shows  the 
nozzle  inlet  conditions  for  the  two  models.  This 
corresponds  to  combustor  exit  flow  in  region  1  obtained  from 
a  cycle  code  using  free  stream  conditions  (6).  Figures  18, 
19,  and  20  prove  the  modified  code  works  for  an 
internal/external  flow  case.  Trends  are  reversed  compared 
to  the  interior  flow  nozzle.  Now  the  imperfect  gas  has  a 
larger  pressure  and  temperature  magnitude  change,  shown  in 
Figures  18  and  19,  respectively.  This  is  due  to  the  values 
of  7  and  Rg^s  (Y=  1.25  and  Rgas=  332.26)  used  for  the  perfect 
gas  model.  This  reversal  demonstrates  the  sensitivity  of 
the  flow  field  solution  to  changes  in  the  specific  heat 
ratio  and  the  gas  composition,  v.iiich  manifests  itself  in 
Rqas.  Though  not  easily  seen  in  Figure  18,  pressure  for  the 
perfect  gas  is  slightly  higher  than  that  for  the  imperfect 
gas.  The  coupling  trend  between  pressure  and  temperature, 
where  temperature  shows  the  larger  diference  betv/een  the 
thermodynamic  models  (Figure  19),  is  the  same  as  that  seen 
for  the  interior  nozzle.  The  modified  code  can  not  yet 
handle  one  mixture  for  the  interior  and  another  for  the 
exterior  flow,  so  no  comparisons  are  made  between  the  two 
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thermodynamic  models  or  with  the  interior  nozzle. 

Central  processing  unit  (CPU)  run  time  was  obtained  for 
the  internal/external  flow  geometry  to  compare  computational 
efficiency.  The  perfect  gas  model  required  2.1  seconds  of 
CPU  run  time  while  the  imperfect  gas  model  finished  the 
solution  in  7.9  seconds.  This  difference  is  a  result  of  the 
iterative  procedures  needed  to  determine  properties  aft  of 
the  waves  for  solution  of  the  Riemann  problem.  Considering 
the  already  remarkable  performance  of  the  first-order 
accurate  FDS  scheme,  this  is  not  a  large  penalty  to  pay  for 
a  more  physically  correct  solution. 

5 . 4  Summary 

Along  with  the  validation  study  of  oblique  shock 
reflections,  the  interior  flow  nozzle,  again,  demonstrates 
the  effect  of  accounting  for  the  caloric  imperfections  of 
high- temperature  gases.  Since  the  net  thrust  margins 
required  to  launch  a  NASP-type  vehicle  are  very  small,  the 
thrust  difference,  found  in  the  proceeding  nozzle  study,  are 
significant.  Realistic  first-order  analysis  of  hypersonic 
nozzles  is  important  to  final  nozzle  design. 

The  imperfect  gas  model  does  not  account  for  all  the 
physics  of  the  flov/.  The  frozen  or  non-reacting  flow 
assumption  made  here  does  not  account  for  dissociation 
losses  nor  energy  recovery  due  to  recombination.  But  a 
study  by  Snelling  (7),  which  implemented  a  finite-rate 
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chemistry  model  into  high-expansion  rate  flow  fields  and 
nozzle  geometries  similar  to  ones  used  here,  concluded  that 
there  is  little  difference  between  the  frozen  flow 
assumption  and  the  finite-rate  chemistry  model  in  solving 
the  flow  field.  This  supports  the  fact  that  the  improved 
thermodynamic  model  can  ir-  're  realistically  model  the  actual 
physics  of  such  flow  fields.  Where  differences  turn  out  to 
be  small  between  thermodynamic  models,  the  first-order 
accurate,  imperfect  gas  FDS  method  should  be  significantly 
more  efficient  than  a  finite-rate  algorithm  in  terms  of  CPU 
time . 
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Figure  10  Parabolic  nozzle  contour 
(3:46) 
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Figure  12  Temperature  along  nozzle 
wall,  IC  1 
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Figure  13  Interior  nozzle  thrust,  IC  1 
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Figure  15  Interior  nozzle  temperature 
IC  2 
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Figure  16  Interior  nozzle  thrust 
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Figure  17  Internal  and  external  nozzle  flow  and 
geometry  (3:60) 
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Figure  18  Pressure  along  nozzle  wall 
for  internal/external  flow  nozzle 
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Figure  19  Temperature  along  nozzle  wall 
for  internal/external  flow  nozzle 
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Figure  20  Integrated  thrust  along 
nozzle  wall  for  internal/external  nozzle 
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Table  2  Initial  conditions  for  internal 
flow  nozzle 


IC  1 

IC  2 

flow  angle  (9) 

0° 

0° 

Mach  number 

6.0 

3.465 

static  pressure  (N/m^) 

1696.4 

200703 

static  temperature  (K) 

273.23 

3079.1 

Table  3  Exterior  initial  conditions  for  the 
internal/external  flow  nozzle 


Perfect  Gas 

Imperfect  Gas 

Mach  number 

8.16 

8.42 

static  pressure  (N/m^) 

4095.6 

4068 . 9 

static  temperature  (K) 

800.2 

776.4 

specific  heat  ratio,  y 

1.4 

1.4 

gas  constant,  Rgas  (J/kg/K) 

287 

287 
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Table  4  Interior  initial  conditions  for  the 
internal/external  flow  nozzle 


Perfect  Gas 

Imperfect  Gas 

Mach  number 

3.465 

3.465 

static  pressure  (N/m^) 

200703 

200703 

static  temperature  (K) 

3079.1 

3079 . 1 

specific  heat  ratio,  y 

1.25 

1.28 

gas  constant,  (J/kg/K) 

332.26 

287.0 
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VI .  Summary  and  Recommendations 


6 . 1  Summary 

The  thermodynamic  model  of  a  previously  existing  CFD 
algorithm  has  been  improyed  by  implementation  of  a  thermally 
perfect  and  calorically  imperfect  gas  model.  By  modifying 
the  thermodynamic  model  and  maintaining  the  desirable 
characteristics  of  the  original  marching  scheme,  a  more 
physically  correct  solution  of  supersonic  planar  nozzle 
flow  is  ayaiJable. 

The  imperfect  gas  model  was  shown  to  be  nearly 
identical  to  exact  solutions  that  were  used  in  the  oblique 
shock  reflection  yalidation  studies.  Physically  real  flow 
geometries,  as  well  as  the  yalidation  studies,  reyealed  that 
the  perfect  gas  model  oyerpredicted  magnitude  changes  in 
performance  quantities  as  pressure  and  temperature  were 
increased.  CPU  run  times  were  3.75  times  longer  for  the 
imperfect  gas  but  the  adyantages  of  a  more  physically 
correct  solution  outweigh  the  reduction  in  computational 
efficiency.  .The  full  potential  of  the  imperfect  gas  model 
will  be  realized  when  combustion  species  can  be  run  for 
internal  nozzle  flow  and  air  for  external  flow. 

In  the  design  and  optimization  of  a  hypersonic  nozzle, 
the  differences  between  thermodynamic  models  will  show  up  in 
the  size  of  the  nozzle  required  to  proyide  a  specified 
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thrust.  In  the  design  of  a  NASP-type  vehicle,  these  results 
will,  in  turn,  influence  engine  size,  fuel  onboard,  vehicle 
size,  and  so  on.  Therefore,  new  and  more  accurate 
thermodynamic  models  will  be  crucial  to  optimizing  the 
entire  vehicle  design. 

6 . 2  Recommendations 

Inevitably,  follow-on  research  topics  arise  as  a  by¬ 
product  of  such  a  study.  Some  future  areas  of  interest  to 
accompany  this  research  would  be; 

1.  Further  modification  of  the  thermodynamic 
model  to  include  finite  rate  chemistry  and 
equilibrium  and  nonequilibrium  flows. 

2.  Experimental  analysis  of  a  simple  nozzle 
geometry  to  validate  the  imperfect  gas  model  in 
the  FDS  algorithm. 

3.  An  accuracy  and  efficiency  comparison  i  '  the 
modified  FDS  algorithm  with  a  parabolized 
Navier- Stokes  code. 

4.  Design,  optimization,  and  comparison  of  planar 
supersonic  nozzles  using  the  perfect  and  imperfect 
thermodynamic  models. 
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Appendix  A:  Thermodynamic  Model 


A. 1  Introduction 

The  complexity  of  a  thermodynamic  model  depends  on  how 
much  physics  of  the  flow  one  wishes  to  account  for.  A 
perfect  gas  is  the  simplest,  where  specific  heats,  c„  and 
Cp,  and  hence,  the  specific  heat  ratio,  y,  remain  constant 
throughout  the  flow.  In  this  case,  enthalpy,  h,  and 
specific  internal  energy,  6,  are  explicit  functions  of 
temperature  (1:388) : 


h  =  CpT 

(47) 

Q  *  C^T 

(48) 

Again,  the  thermal  equation  of  state  holds  and  for  frozen 
flow,  the  specific  gas  constant,  Rg^s,  remains  constant  since 
the  molecular  weight  of  the  gas  is  fixed: 

n 

n  -  ^universal  fAQ\ 

- mT*  '  ’ 

For  the  perfect  gas  assumption,  many  basic  thermodynamic 
relations  can  be  solved  for  explicitly  and  shock  wave 
analysis  is  comparatively  straightforward. 

At  the  other  end  of  the  spectrum  is  a  nonequilibrium, 
chemically  reacting  gas.  This  model  accounts  for  most  high- 

temperature  gas  effects.  Flow  properties  are  not  only  a 

function  of  temperature  but  also  pressure  and  time.  Time  is 
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an  independent  variable  since  reactions  have  not  reached 
steady  state.  The  effect  of  pressure  is  to  inhibit 
reactions  at  higher  pressures  and  promote  them  at  lower 
pressures  ( 1 : 374 ) . 

Dissociation  is  another  high- temperature  effect. 
Dissociation,  which  alters  the  molecular  weight  of  the  gas, 
occurs  at  around  2500  K  for  Oj,  and  begins  at  4000  K  for  Nj- 
Van  der  Waals'  effect  also  comes  into  play  as  intermolecular 
forces  become  important  (1:390). 

A . 1  The  Imperfect  Gas  Model 

When  categorized,  by  assumptions,  among  the  various 
thermodynamic  models,  an  imperfect  gas  assumes  that 
intermolecular  forces  are  negligible.  With  respect  to  a 
perfect  gas,  an  imperfect  gas  is  the  next  step  in 
complexity.  The  difference  between  these  two  models  is  a 
function  of  which  mode  of  internal  energy  is  excited.  When 
dealing  with  monatomic  gases,  specific  heats  are  basically 
independent  of  temperature  since  only  the  translational 
energy  mode  is  available  (neglecting  electronic  energy) . 
Since  most  flow  field  gases  are  diatomic,  the  translational, 
vibrational,  and  rotational  energy  modes  are  available,  even 
at  room  temperature  (5:222).  For  air  as  a  perfect  gas,  it 
is  assumed  that  internal  energy  is  only  a  function  of 
translational  and  rotational  degrees  of  freedom  (in  this 
discussion  of  energy  modes,  e  will  represent  specific 
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internal  energy  instead  of  fi) : 

^perfect  ~  ^cr  *  ^zoc 

=  —RT  +  RT 
2 


(50) 


Recalling  that  e=CvT,  the  specific  heat  at  constant  volume 
yields  Cv=5/2R.  For  a  thermally  perfect  gas: 

Cp  =  R 


(52) 


Where  the  familiar  quantity,  7=  1.4,  is  obtained  for  air  as 
a  perfect  gas  (9:124,133). 

When  temperature  reaches  600  K,  vibrational  energy  is 
no  longer  negligible  (1:441)  and  a  new  internal  energy  mode 
must  be  added.  From  a  statistical  mechanics  approach,  this 
component  is  represented,  without  further  proof,  as  (1:439): 


hv/kT  „„ 

g/jv/fcT_2 


(53) 


where  h  is  Planck's  constant,  k  is  Boltzmann's  constant,  and 
V  is  the  fundamental  vibrational  frequency  of  the  molecule. 
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The  imperfect  gas  model  accounts  for  this  vibrational  mode. 

In  summary,  four  modes  of  energy  each  contribute  a 
portion  to  the  total  internal  energy  (1:440): 

^  I  *  %1  (5«) 

Comparing  this  equation  to  eq.  (50)a,  it  is  obvious  that  the 
energy  of  gaseous  molecules  is  spread  over  more  modes  for  an 
imperfect  gas  (three  modes)  than  for  a  perfect  gas  (two 
modes) .  For  perfect  gas  flow  across  a  shock,  the  kinetic 
energy  in  front  of  the  shock  is  converted  to  translational 
and  rotational  energy  aft  of  the  shock.  For  imperfect  gas 
flow,  the  vibrational  energy  mode  can  accommodate  some  of 
this  energy  conversion  across  the  shock,  lessening  the 
amount  of  energy  going  into  the  translational  and  rotational 
modes,  relative  to  a  perfect  gas.  This  effect  is  important 
at  temperatures  above  600  K.  Since  temperature  is  a  direct 
measurement  of  kinetic  (translational)  energy,  a  perfect  gas 
tends  to  overpredict  temperature  (1:510). 

A . 2  Implementing  Imperfect  Gas  Assumptions 

Determining  where  modifications  to  Doty's  algorithm  (3) 
are  required  is  a  matter  of  locating  what  calculations  are 
based  on  a  perfect  gas  model.  The  most  obvious  changes  are 
in  the  solution  of  the  Riemann  problem.  These  are  basically 
compression  or  expansion  wave  problems  and  applying 
imperfect  gas  assumptions  is  v/ell  documented  in  text 
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materials  (10).  For  the  linear-approximate  solution  that  is 
being  used  for  this  study,  the  method  for  finding  pressure 
aft  of  the  discontinuity  remains  the  same.  What  is  less 
straightforv/ard  is  how  to  calculate  properties,  such  as 
entropy  and  enthalpy,  in  the  aft  shock  region  based  on 
variable  specific  heats.  Since  decoding  of  the  E  and  F 
vectors  is  based  on  constant  y,  this  must  also  change. 

As  is  mentioned  in  the  beginning,  a  desirable 
capability  is  to  handle  a  multiple  species  gas.  The  current 
perfect  gas  code  simulates  air  (physically  based  on  a 
combination  of  nitrogen,  oxygen,  and  argon)  as  the  working 
fluid.  The  only  way  to  alter  the  composition  of  the  gas  is 
through  the  gas  constant,  where  R^n^ymole-wt .  It 

does  not  have  the  capability  to  handle  individual  gas 
species.  To  allow  greater  flexibility,  a  multiple  species 
gas  would  allow  this  code  to  handle  the  output  from  a 
separate  cycle  code  analysis  program.  The  combustion 
products  of  a  SCRAMJET  engine  will  consist  of  mostly  water 
and  nitrogen  which  are  quite  different  from  air  as  far  as 
molecular  weight  goes,  and  hence,  will  behave  differently. 
Doty's  code  can  handle  a  nonuniform  input  from  a  combustor 
exit  and  this  capability  would  be  enhanced  with  true 
combustor  products  rather  than  air. 

A . 3  Complex  Chemical  Mixtures 

A  FORTRAN  code,  referred  to  here  as  CET89,  by  Gordon 
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and  McBride  (4),  outputs,  among  other  things,  the  least 
squares  coefficients  of  specific  heats,  enthalpy,  and  an 
entropy  parameter  as  functions  of  temperature.  Data  from 
JANAF  thermochemical  data  has  been  curve  fitted  for 
approximately  400  species  of  gas  using  the  method  of  least 
squares  to  the  obtain  coefficients.  There  are  two  sets  of 
coefficients  for  each  species  that  correspond  to  two 
temperature  ranges,  300-1000  K  and  1000-5000  K.  The 
calculation  of  the  three  properties  are  listed  below. 

A. 3.1  Specific  Heat,  Cp. 

Specific  heat  for  a  mixture  of  gases  is  represented  by 
the  following  summation  (10:50): 

where  n  is  a  particular  gas  species  and  C,  is  the  mass 
fraction  of  the  species.  In  terms  of  the  least  squares 
coefficients  Cp  is  written,  on  a  mass  basis,  as  (4:32): 

Cp^  =  (a  +  Jbr  +  cT  2  +  dr-  +  (56) 


A. 3. 2  Enthalpy,  h. 

Enthalpy  for  a  thermally  perfect  gas  is  given  by  the 
following  expression  (10:50): 
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(57) 


h  =  fc  dt 

Substituting  eq.  (56)  into  this  equation  yields  the  least 
squares  coefficients  form  of  enthalpy,  for  the  i'"’  species, 
on  a  mass  basis  (4:32): 


^  ij*  2 
2 


d  rp  4 

4 


h 

“298  1  ^ 

MW  ^ 


where  hjjg  is  the  mass  basis  value  of  enthalpy  at  T=  298  K. 


A. 3. 3  Entropy, _ s. 

Entropy,  for  a  thermally  perfect  gas,  is  given  by  the 
following  expression  (10:51): 


=  /' 

•>  Co 


=  <})  -  Rln- 


(59) 


Again,  substituting  eq ,  (56)  into  the  integral  portion  of 
this  equation  yields  the  least  squares  coefficient  form  of 
the  entropy  parameter,  <{), ,  for  the  i'^’’  species  (4:32): 

4.,.  =  (4>„  +  ainr  +  bT  *  (60) 

Appendix  B  and  Appendix  D  show  how  these  expressions 
are  used  to  implement  the  imperfect  gas  model. 

A . 4  Code  Modifications 

The  original  code  consisted  of  46  subroutines.  Three 


57 


ware  added  and  eight  were  modified.  The  following  is  a 
summary  of  the  altered  subroutines.  Appendix  E  contains  the 
three  additional  subroutines  and  imperfect  gas  oblique  shock 
wave  solver. 

Subroutine  INITEM:  initializes  the  arrays  and 
variables.  Values  of  coefficients  from  CET89  output 
are  listed  here. 

Subroutine  INPUT:  read  statements  for  input  deck. 

Mass  fractions  are  read  in  here. 

Subroutine  INVALU:  initial  value  line  for  inlet  to 
nozzle  used  to  begin  Riemann  problem.  Speed  of  sound 
and  density  require  Y(T)  and  gas  constant,  Rga^,  for 
mixture . 

Subroutine  OUTINP:  outputs  the  input  deck. 

Subroutine  RMCONT:  determines  which  Riemann  solver  to 
use.  Updated  z  coefficient  requires  Y(T) . 

Subroutine  RMAPAP:  solution  of  the  Riemann  problem. 

Subroutine  UPWINDl:  marches  the  solution  in  the  x 
direction . 

Subroutine  BNDWAV:  marches  the  solution  in  x  direction 
along  boundaries. 
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Appendix  B:  Solution  of  the  Riemann  Problem 


This  discussion  of  the  linearized-approximate  method 
applied  to  the  solution  of  a  perfect  gas  Riemann  problem  is 
summarized  from  Doty's  research  (3).  The  reader  is  urged  to 
consult  reference  (3)  for  more  complete  details. 

As  previously  mentioned,  the  linearized-approximate 
method  treats  waves  as  isentropic  expansions  and 
compressions.  By  linearizing  the  Prandtl-Meyer  equations,  a 
closed-form  equation  results  that  requires  no  iteration. 
Figure  21  illustrates  the  setup  of  the  Riemann  proble.m. 

B . 1  Linearized-Approximate ,  Perfect  Gas 

For  isentropic,  planar,  steady  flow,  the  differential 
form  of  the  compatibility  relations,  valid  along  Mach  lines, 
is : 

7^2-1  dP  ±  pV^dQ  =  0  ) 

where  the  velocity  magnitude  and  flow  angle  are: 

y2  =  jj2  ^^2  (  62  ) 

0  =  tan'^  ( v/u)  ( ) 

Through  substitution,  manipulation,  and  realizing  that 
p=7P/a^,  eq.  (61)  is  recast  as: 

.3P  t  d(v,/u)  .  0  (64) 
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Two  definitions  and  a  relationship  are  introduced: 


(65) 


o  =  v/u  (66  ) 

d[ln(P)]  =  -^  (67) 

Eq.  (61)  is  then  written  more  compactly  as: 

d[ln(P)]  ±  (z)da  =  0  (68) 

Linearizing  this  yields: 

A[ln(P)]  ±  (z)Ao  =  0  (69) 

Referring  to  Figure  21,  a  positive  wave,  wave  (3),  is 
required  to  pass  information  from  region  (0)  to  region  (2). 
Using  the  positive  sign  from  eq.  (69)  yields: 

( [ln(P)  ]  2“  [ln(P)  Iq)  +  (Zq)  (Oj-Oq)  =0  (70) 

Rearranging  eq.  (70)  yields: 

[ln(P)  ]  2+ (Zq^  =  [ln(P)  ]„+ (Zq)  Oq  (71) 

A  similar  process  is  performed  for  regions  (6)  and  (4): 

[ln(P)  ]4+ (Zg)  O4  =  [In  (P)  )g+ (Zg)  Og  (72) 

Recall  that  the  pressure  and  flow  slope  must  match  across 
the  contact  surface: 
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(73) 


=  O2 

P4  =  P2  (74) 

Eqs .  (71),  (72),  (73),  and  (74)  represent  four  equations 
and  four  unknowns  and  may  be  solved  in  closed  form. 
Substituting  eqs.  (73)  and  (74)  into  eq .  (71)  and 
rearranging  gives: 

[ln(P)  ]  4+  (Zg)  O4  =  [In  (P)  ]  g+  (Zg)  Og  (75) 

Solving  for  O4  yields: 

_  _  [ln(P)  ]g-[ln(P)]g+(z6)06+(Zg)0g 

Using  this  value  in  the  solution  of  P4,  in  eq.  (72),  gives: 

P4  =  exp(  [ln(P)  ]6+Zg(04-06) )  (77) 

Isentropic  relation.,  can  be  used  to  obtain  density  and  speed 
of  sound  in  regions  (2)  and  (4)  (only  region  (2)  is  shown 
here)  : 

Pi  = 

Po 

32  =  (YP2/P2]'/'  (79) 

Conservation  of  stagnation  enthalpy  across  a  wave  is 
used  to  obtain  velocity  (recall  that  h=CpT)  : 


l/Y 


(78) 
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(80) 


_ 

^4 

1+ 

7; 

Eq.  (80)  can  be  solved  explicitly  for  M4 ;  and  the  velocity 
components,  in  terms  of  flow  slope,  O4,  are: 


U4  =  Af^a^cos  [tan'Mo^)  ] 

(81) 

V4  =  Af^a^sinCtan'Mo^)  ] 

(82) 

A  similar  process  is  done  for  region  (2). 

B . 2  Linearized-approximate ,  Imperfect  Gas 

Pressure  in  region  (4),  P4,  from  eq.  (77),  is  also 
valid  for  an  imperfect  gas.  Since  pressure  is  the  only 
property  known  in  region  (4),  a  relation  is  needed  that  does 
not  rely  on  a  calorically  perfect  form  of  conservation  of 
stagnation  enthalpy. 

To  determine  temperature,  relations  are  used  that  take 
advantage  of  the  isentropic  nature  of  the  flow.  For  a 
thermally  perfect  gas,  the  differential  quantity  of  entropy 
in  a  system  is  (10:45): 

ds  =  Cp^  -  (83) 

Integrating  this  yields  (10:58): 

s  =  —  -  Rln—  =  <f>  -  Rln—  (84) 

•'Co  -P  C  Po 

where  the  entropy  parameter,  (j),  is  defined  as  (10:58): 
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(85) 


Substituting  eq.  (56),  on  a  molar  basis,  into  eq.  (85) 
yields  the  already  integrated  equation  (10:58): 


4) 


4)  +  alnt  +  jbt  + 


ct' 


dt- 


et^ 


R 


(86) 


Appendix  A  shows  how  the  entropy  parameter,  (|),  can  be 
determined.  A  change  in  entropy  is  found  from  eq.  (84)  and 
written  for  regions  (4)  and  (6)  as  (10:58): 

^4  -  ^6  =  4>4  -  4>6  -  In  j^-^j  (87) 


For  isentropic  flow. 


As=0  and  eq.  (87) 


*  R  In 


is  rewritten  as: 


(88) 


where  ({)g  is  found  using  eq.  (86)  in  terms  of  the  known 
temperature  in  region  (6).  The  entropy  parameter  in  region 
(4),  <{)4,  is  the  only  unknown  is  this  equation.  Using  eq. 
(86),  temperature  in  region  (4),  Tj,  is  solved  for 
iteratively  using  the  secant  method.  Subroutine  ENTROPY 
performs  this  operation  and  can  be  found  in  Appendix  E. 

With  temperature  and  pressure  in  region  (4),  density  can  be 
determined  using  the  thermal  equation  of  state. 

As  in  the  perfect  gas  case,  conservation  of  stagnation 
enthalpy  is  used  to  find  the  velocity  components,  but  in  a 
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different  form.  Stagnation  enthalpy  is  determined  in  region 
(4)  as; 


=  hg  +  +  (89) 

where  static  enthalpy,  hg,  is  a  function  of  Tg  and  is  found 
using  the  methods  described  in  Appendix  A.  Across  wave  (3), 
stagnation  enthalpy  in  region  (4)  is  found  using  the 
conservation  of  energy  relationship,  h^g  =  h^ ,  Subroutine 
ENTHGAM  performs  the  calculation  of  enthalpy  and  ratio  of 
specific  heats  based  on  temperature  and  can  be  found  in 
Appendix  L.  Recalling  thut  o=v/u,  stagnation  eiiLiiaipy  can 
be  written  as: 


(I+O4) 

Solving  for  u^  gives: 


(90) 


(91) 


The  y  component  of  velocity  is  then  computed  as: 

.V4  =  U4O4  (92) 

Mach  number  can  now  be  computed: 
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Appendix  C:  Numerical  Algorithm 
C . 1  Introduction 

The  methodology  of  flux-difference-splitting  for 
supersonic  planar  flow  and  the  first-order  accurate  marching 
scheme  is  presented  in  this  appendix.  This  is  a  summary  of 
the  numerical  schemes  and  the  reader  is  referred  to 
reference  (3)  for  further  details. 

C . 2  Riemann  Fluxes 

After  solving  the  Riemann  problem  by  one  of  the  three 
.methods  described,  the  Riemann  fluxes  are  computed. 

Figure  22  shows  a  schematic  of  the  flux  differencing  and 
splitting;  The  solution  of  the  Riemann  problem  in  Appendix  B 
is  the  basis  for  computing  the  fluxes. 

To  illustrate,  consider  the  flux  vector,  El,  for  the 


Riemann  regions  0,  2,  4,  6: 

{E1)q  =  PqUq  (94) 

{£1)2  =  P2U2  (95) 

iEl)^  =  P4U4  (96) 

(El  )g  =  pgUg  (97) 


The  Riemann  flux-differences  are  the  differences  of  the  El 
vector  taken  across  waves  (1),  (2),  and  (3): 
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(98) 

=  {E1),-(E1)^ 

(99) 

{dED,^^^,  =  (E1)^-{E1)^ 

(100) 

The  sum  of  the  contributions  across  all  three  waves  gives 
the  total  contribution  at  Riemann  node  j+1/2  (3:184): 

(dEl'^.i/2  =  (101) 

Riemann  node  j-1/2  is  handled  in  a  similar  fashion,  as  are 
the  other  E  and  F  vector  components. 

C . 3  Splitting  the  Flux  Differences 

With  the  Riemann  problem  solved  and  the  flux 
differences  taken,  the  differences  are  now  split  to 
determine  which  information  from  Riemann  nodes  j+1/2  and  j- 
1/2  will  reach  node  j  at  the  next  plane,  i+1  (3*183). 

Figure  22  illustrates  the  flux  difference  and  splitting. 

Streamlines  and  characteristics  determine  which 
direction  to  send  the  differenced  fluxes.  As  described  by 
Doty  from  Chakravarthy ' s  report,  Godunov  and  Osher  developed 
identical  methods  for  splitting  the  flux  differences,  in 
unsteady  flow,  based  on  the  sign  of  the  slopes  of  the  waves 
emanating  from  the  Riemann  node  (3:184).  This  concept  was 
extended  to  steady,  two-dimensional  flow.  For  example,  at 
node  j-1/2,  if  the  slope  of  wave.s  (1),  (2),  or  (3)  is 
positive,  the  flux  difference  cf  the  E  and  F  vectors  across 
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that  wave(s)  contributes  to  the  solution  at  node  i+l,] 
(3:185).  A  similar  process  occurs  at  Riemann  node  j+1/2 
where  this  time  only  flux  differences  with  a  negative  slope 
contribute  to  the  solution  at  i+l,j.  In  equation  form, 
summing  the  positive  and  negative  contributions  at  the 
Riemann  nodes  j+1/2  and  j-1/2,  respectively,  yields  (3:186): 

dEj.,„  =  -  KT?!)'  *  ' 

dE;.,,,  =  [dE“r;iT  •  [d£A7ii‘  *  [d£7i™T 

where  the  positive  or  negative  signs  denote  positively  or 
negatively  sloped  flux  differences.  All,  none,  or  some  of 
the  waves  can  contribute  to  the  solution,  depending  on  their 
slopes.  The  same  procedure  is  applied  to  the  F  vector.  For 
further  details,  refer  to  Appendix  J  of  reference  (3). 


C . 4  First-Order  Accurate  Marching  Scheme 

The  governing  equations  for  supersonic  flov;  are 
hyperbolic  and  can  thus  be  solved  by  a  marching  technique. 
The  transformed  governing  equations  are  developed  in 
Appendix  D  of  reference  (3)  and  are  repeated  here  without 
proof  ( 3 : 188  )  : 


d{E)  _  d{E)  d{F) 

ac  '  ^  ari  dr\ 


(104) 


In  order  to  march  the  flov;  in  the  x  or  C  direction  a 
combination  of  finite  difference  and  flux  difference 
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approximations  are  used  for  the  partial  derivatives  of  E  and 
F  in  eq.  (104).  Substituting  these  approximations  directly 
into  eq.  (104)  yields: 


A^{E)  AAB)  AAF) 

-at 


(105) 


where  A^  is  the  finite  difference  operator  for  the  C 


direction  and  A,  is  the  flux  difference  operator  for  the  ri 
direction  (3:188).  For  convenience,  At]  is  chosen  to  be 


unity.  Rearranging  yields: 


A^iE)  =  -ACn^Aj(E)  -ACT)yAj.(F)  (106) 

The  solution  is  advanced  in  the  x  direction  using  a 
finite  difference  operator  for  A,(E)  in  eq.  (106): 


Aj;{E)  =  -  e}  (107) 

Substituting  this  into  eq.  (106)  and  rearranging  yields 
(3 : 189) : 


=  E^  -  AC^^AJ{E)  -  ACTiyAj(F)  (108) 

As  discussed  earlier,  the  solution  uses  negatively 
biased  information  from  Rremann  node  j+1/2  and  positively 
biased  information  from  Riemann  node  j-1/2.  Recalling  the 
results  of  eqs .  (102)  and  (103),  the  first-order  accurate 
FDS  representations  for  A^E  and  A^F  operators  are  (3:190): 


Aj(E)  -  {ciEj.^/2 


(109) 
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Aj.(F)  =  (dFj.i/j  -  dFj'.i/a)  (110) 

Substituting  eqs .  (109)  and  (110)  into  eq.  (108)  yields  the 
first-order  accurate  FDS  approximation  to  the  solution  of 
the  governing  equations  (3:190): 

(111) 

=  e}  -  A(Tl4dF;.i/2  +  dFj^i/a)  -  ACtly[dF;.a/2  +  dF^^^/a) 
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Appendix  D;  Decoding  the  Solution 


D . 1  Introduction 

After  marching  the  solution  to  the  next  plane,  i+1,  the 
newly  calculated  E  and  F  vectors  must  be  decoded  into 
primitive  variables  so  the  Riemann  problem  may  be  solved. 
This  will  then  complete  the  "loop"  needed  to  solve  the  flow 
field  from  one  plane  to  next. 

D . 2  Decoding,  Perfect  Gas 

Decoding  the  E  vector  requires  that  all  components  be 
solved  simultaneously.  The  E  vector  with  its  four 
components  are  known  from  marching  to  the  i+1  plane.  It  is 
repeated  here  for  convenience: 

pu 

E=  (112) 

puv 

u  (pe+P) 

The  pe  term  in  the  E4  component  can  be  written  as: 

pe  =  pu  +  -ip(u^  +  v2)  (113) 

For  a  perfect  gas,  specific  internal  energy,  6,  can  be 
written  as: 

a  =  c^T  (114) 

The  following  relation  holds  for  a  thermally  perfect  gas 


72 


V 


(10:44)  : 


,  _  ^gas 

'''  ~  Y-1 


(115) 


Recalling  the  thermal  equation  of  state,  substituting  eq. 
(115)  into  eq.  (114),  and  substituting  that  result  into  eq. 
(113)  yields: 


^  +  |p(u2+v2) 


(116) 


This  is  the  fifth  equation  to  close  the  system  of  five 
unknowns.  Eq.  (116)  can  be  substituted  into  the  E4 
component  to  give: 


E4  =  u  [  (— ^+  — p  (u^  +  v^) ) +P] 
y-l  2 


(117) 


Expanding  and  rearranging  eq.  (117)  gives: 


E4  -  UP  — ^1+  — pu (u^  + v^) 
Y-lJ  2 


(118) 


Recalling  that  El=  pu,  E2,  E3,  and  E4  can  be  written  as: 


E2  =  {ED  u+P 


(119) 


E3  =  (EDv 


(120) 


E4  =  u  +  l  (El)  (uD  (El)  ( v2) 

Y-lJ  2  2 


(121) 


Solving  eq.  (120)  for  the  y  component  of  velocity,  v,  gives 
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u 


El 

El 


(122) 


Eq.  (119)  can  be  solved  for  P: 

P  =  E2- {El)  u  ( 123  ) 

Eqs .  (122)  and  (123)  can  be  substituted  into  eq.  (121)  to 

yield : 

E4  =  u\-l—][{E2) -{ED  u] +- (El)  {u^)  +-  {E1)\—Y  (124) 

y-lj  2  2  [El 

These  manipulations  have  resulted  in  an  E4  component  that  is 
a  function  of  the  known  E  vector  components  and  the  one 
unknown,  u,  and  can  now  be  cast  as  a  quadratic  equation: 

— (PI)]  -f— I- (£)2)]  u  +  {E4)-E.1E111  =0  (125) 

[2(y-1)  “  [y-1  ]  [  2  {ED 

au^+bu+c  =  o  (126) 

This  quadratic  equation  can  now  be  used  to  find  the  x 
component  of  velocity,  u,  where  it  has  been  determined  that 
the  positive  value  is  the  correct  root  (3:222).  The 
primitive  variables  are: 
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-Jb  +  \/b^  -  4ac 

2a 

(127) 

V  =  ^ 

El 

(128) 

P  =  E2-{E1)U 

(129) 

p  , 

u 

(130) 

pe  =  (D^  +  v^) 

(131) 

This  completes  the  decoding  process  for  a  perfect  gas. 

D . 3  Decoding,  Imperfect  Gas 

Decoding  the  E  vector  for  an  imperfect  gas  requires 
that  the  pe  term  of  the  E4  component  be  consistent  with  the 
imperfect  gas  assumption.  Recall  eq.  (113)  ,  rewritten  here 
for  convenience: 

pe  =  pu+-ip(u^  +  v^)  (132) 

Specific  internal  energy,  dl,  for  an  imperfect  gas  is  given 
as  (10:58): 

G  ^  h  -  RT  ( 133 ) 

Substituting  eq.  (133)  into  eq.  (132)  yields  the  thermally 
perfect  form  of  total  internal  energy: 


75 


(134) 


pe  =  ph  -  P  +  p-|  (u^  +  v^) 


Substituting  eq.  (134)  into  the  E4  component  and  canceling 
the  pressure  terms  gives: 


E4  =  u[(ph  -  P  +  —  p(u^  +  v^)  +P] 


(135) 


E4  =  puh  +  -ipu(u^  +  v2)  (136) 

Substituting  the  El  component  into  eq.  (136)  yields  the 
final  form  of  the  E4  component: 

E4  =  h{El)  +  -  {ED  +  iSMill  (137) 

2  2  El 

This  equation  is  now  a  function  of  the  x  component  of 
velocity,  u,  and  enthalpy,  which  itself  is  a  function  of 
temperature  from  the  least  squares  coefficient.  An 
expression  for  the  velocity  component,  u,  must  be  found  that 
IS  a  function  of  temperature  so  as  to  obtain  an  equation  for 
E4  that  will,  subsequently,  also  have  one  unknown, 
temperature.  Start  by  solving  the  El  component  for  u:  u= 
El/p.  Substitute  the  thermal  equation  of  state  for  p: 

u  =  (138) 

P 

Solve  the  E2  component  for  P  and  substitute  the  result  into 
eq .  ( 138 ) : 
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u 


[EDRT 
E2-[E1)  u 


(139) 


Rearrange  eq.  (139)  into  a  quadratic  equation  in  the  x 
component  of  velocity,  u: 

(El)  +  {E2)U  +  (El)RT  =  0  (140) 


U 


2  (El) 


(141) 


From  this  manipulation,  it  is  seen  that  u  is  now  a  function 
of  temperature.  As  was  desired,  the  E4  component  is  now 
also  a  function  of  only  temperature.  Temperature  can  now  be 
found  by  an  iterative  technique,  such  as  the  secant  method. 
With  temperature  determined,  u  can  be  obtained  by  the 
quadratic  equation.  Similar  to  the  perfect  gas  case, 
pressure  can  be  obtained  from  eq.  (129)  and  density  is 
obtained  from  eq.  (130). 

The  decoding  process  for  the  imperfect  gas  model  is  now 
complete  and  accounts  the  variable  specific  heats  as  a 
function  of  temperature. 
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Appendix  E;  New  Subroutines 


This  appendix  contains  the  three  additional  subroutines 
ENTROPY,  ENTHGAM,  and  IMPERFECT,  and  the  code  for  solving 
the  imperfect  gas  oblique  shock  wave  problem.  Modifications 
to  the  complete  code  are  denoted  by  *  in  FORTRAN  column  one. 
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rHCNfO-cJ'LO'.Or'OO 


subroutine  entropy  (pl,p2,tt,t) 


c 

c  computes  temperature  after  an  isentropic  wave  based 
c  on  the  pressure  in  the  regions  before  and  after  wave 

c  and  temperature  in  the  region  prior  to  the  wave, 

c 


common/coefficient/ 

1  alN2,blN2,clN2,dlN2,elN2,plN2, hlN2, 

2  ahN2,bhN2,chN2/dhN2,ehN2,phN2, hhN2, 

3  al02,bl02,cl02,dl02,el02,pl02, hl02, 

4  ah02 , bh02 , ch02 , dh02 , eh02 , ph02 , hh02 , 

5  alAr , plAr , hlAr , 

6  ahAr,phAr,hhAr/ 

7  alH,blH,clH,dlH,elH,plH,hlH, 

8  ahH, bhH, chH, dhH , ehH, phH, hhH, 

9  alOH,blOH,clOH,dlOH,elOH,plOH,hlOH, 

9  ahOH , bhOH , chOH , dhOH , ehOH , phOH , hhOH , 

1  alH2,blH2,clH2,dlH2,elH2,plH2,hlH2, 

2  ahH2 , bhH2 , chH2 , dhH2 , ehH2 , phH2 , hhH2 , 

3  alNO,blNO,clNO,dlNO,elNO,plNO,hlNO, 

4  ahNO, bhNO, chNO,  dhNO, ehNO, phNO, hhNO, 

5  alO,blO,clO,dlO,elO,plO, hlO, 

6  ahO,bhO,chO,dhO,ehO,phO, hhO, 

7  alH20,blH20,  clH20,dlH20, elH20, plH20, hlH20, 

8  ahH20, bhH20, chH20, dhH20, ehH20, phH20, hhH20, 

9  lN2,102,lAr,lH,iOH,lH2,lMO,10, 1H20 

common/species/ 

xmsf  N2  ,  xmwN2 ,  xmsf02  ,  xmv/02  ,  xmsf  Ar ,  xmv/Ar , 
xmsf  H,  xmwH,  xmsfOH ,  xmv/OH,  xmsf  H2  ,  xmwH2  ,  xmsf  NO,  xmv/NO, 
xmsfO,xmwO, xmsf H20, xmwH20, 
runiv, rmix 


iter=0 
itmax=50 
rmtol= . 000001 
tl=tt 

if ( tl . le . 1000 . )  then 
phimixl= 

lN2*(plN2+alN2*log(  tl)-f-blN2*tl  +  clN2*tl*tl/2.  +dlN2* 
tl**3  ./3.+elN2*tl**4./4.  )  *runiv*xmsfN2/xmwN2  -r 
102* (pl02+al02*log( tl ) +bl02*tl+cl02*tl*tl/2 . +dl02* 
tl**3 ./3 . +el02*tl**4 ./4 . ) *runiv*xmsf02/xmw02  + 
lAr* (plAr+alAr*log( tl ) ) *runiv*xmsf Ar/xmwAr+ 

IH* (plH+alH*log ( tl ) ) *runiv*xmsf H/xmwH+ 
10H*(pl0H+al0H*log(tl)  ^-blOH*tl  +  clOH*tl*tl/2  .+diOH* 
tl**3 . /3 . +elOH*tl**4 . /4 . ) *runiv*xmsfOH/xmwOH  + 
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'^av(-n4^0JtOM'X>vo 


lH2*(plH2+alH2*log(tl)+blH2*tl+clH2*tl*tl/2 .+dlH2* 
tl**3  ./3  .  +elH2*tl**4  .  /4  .  )  *runiv*x[nsf  H2/xniwH2  + 
lN0*(plN0+alN0*log(tl)+blN0*tl+clN0*tl*tl/2 .+dlNO* 
tl**3 ./3 . +elN0*tl**4 ./4 . ) *runiv*xmsf NO/xmwNO  + 

10* (plO+alO*log(tl)+blO*ti+clO*tl*tl/2 . +dlO* 
tl**3 ./3 • +elO*tl**4 ./4 . ) *runiv*xmsfO/xmwO  + 
lH20*(plH20+alH20*log(tl)+blH20*tl+clH20*tl*tl/2 . 
+dlH20*tl**3 ./3 . +elH20*tl**4 ./4 - ) 

*runiv*xinsf  H20/xmwH20 

else 

phimixl= 

1  lN2*(phN2+ahN2*log( tl)+bhN2*tl+chN2*tl*tl/2 .+dhN2* 

2  tl**3  ./3  . +ehN2*tl**4  ./4  .  )  *^'-iiiiv*xmsfN2/xmwN2  + 

3  102*(ph02+ah02*log(tl)+bh02*tl+ch02*tl*tl/2 .+dh02* 

4  tl**3  . /3  . +eh02*tl**4  . /4  .  )  *^'Jniv*xmsf02/xmv/02  + 

5  lAr* (phAr+ahAr*log( tl) ) *runiv*xmsfAr/xmwAr  + 

6  lH*(phH+ahH*log(tl) )*runiv*xmsfH/xmwH  + 

7  10H*(ph0H+ah0H*log(tl)+bh0H*tl+ch0H*tl*tl/2.+dh0H* 

8  tl**3  ./3  - +ehOH*tl**4  ./4  .  )  *i^'J>^iv*xmsfOH/xmwOH  + 

9  1H2*  (phH2+ahH2*log(tl)4-bhH2*tl+chH2*tl*tl/2  .  +dhH2* 

9  tl**3  ./3  . +ehH2*tl**4  ./4  .  )  *i'i^'^iv*xmsf  H2/xmwH2  + 

1  INO* (?hN0+ahN0*log(tl)+bhN0*tl+chN0*tl*tl/2 . +dhNO* 

2  tl**3 ./3 • +ehN0*tl**4 ./4 - ) *runiv*xmsf NO/xmwNO  + 

3  10* ( phO+ahO*log ( tl ) +bhO*tl+chO*tl*tl/2 . +dhO* 

4  tl**3  ./3  . +ehO*tl**4  ./4  .  )  *tuniv*xinsfO/xmwO  + 

5  lH20*(phH20+ahH20*log(tl)+bhH20*tl+chH20*tl*tl/2 . 

6  +dhH20*tl**3 ./3 .+ehH20*tl**4 ./4 . ) 

7  *univ*xnisf  H20/xmwH20 
end  if 

phimix2=phimixl+rinix*log(p2/pl ) 
t2=tlT50. 

do  100  iter=l , itmax 
if ( t2 . le . 1000 . )  then 
phimix2g= 

1  lN2*(plN2+alN2*log(t2)+blN2*t2+clN2*t2*t2/2.+dlN-2* 

2  t2**3 ./3 . +elN2*t2**4 ./4 ■ ) *runiv*xmsfN2/xmwN2  + 

3  102*(pl02+al02*log(t2)-f-bl02*t2+cl02*t2*t2/2  .■i-dl02* 

4  t2**3 ./3 • +el02*t2**4 . /4 . ) *runiv*xmsf02/xmw02  + 

5  lAr*(plAr+alAr*log(t2) )  *runiv*xmsfAr/xmv/Ar  + 

6  lH*(plH+alH*log(t2) )*runiv*xmsfH/xmwH  + 

7  lOH*  ( plOH+alOH*log  ( t2 )  •fblOH*t2+clOH*t2*t2/2  .  +dlOH* 

8  t2**3 . /3 - +elOH*t2**4 . /4 . ) *runiv*xmsfOH/xmwOH  + 

9  lH2*(plH2+alH2*log{  t2)-i-blH2*t2+clH2*t2*t2/2  .  +dl[12* 

9  t2**3 ./3 • +elH2*t2**4 ./4 . ) *runiv*xmsf H2/xmwH2  + 

1  INO* ( plNO+alNO*log ( t2 ) +blNO*t2+clNO*t2* t2/2 . +dlNO* 

2  t2**3  . /3  . +elNO*t2**4  ./4  .  )  NO/xmwNO  + 

3  10* (plO+alO*log ( t2 ) +blO*t2+clO*t2* t2/2 . +dlO* 

4  t2**3 . /3 - +elO*t2**4 ./4 ■ ) *runiv*xmsfO/xmwO  + 
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5  1H20* (plH20+alH20*log(t2)+blH20*t2+clH20*t2*t2/2 . 

6  +dlH20*t2**3 ./3 .+elH20*t2**4 ./4 • ) 

7  *runiv*xmsf H20/xmwH20 
else 

phimix2g= 

1  lN2*(phN2+ahN2*log(t2)+bhN2*t2+chN2*t2*t2/2 .+dhN2* 

2  t2**3  ./3  ■ +ehN2*t2**4  ./4  •  )  *j^iJniv*xmsf  N2/xmwN2  + 

3  102*(ph02+ah02*log(t2)+bh02*t2+ch02*t2*t2/2 . +dh02* 

4  t2**3 ./3 • +eh02*t2**4 ./4 ■ ) *runiv*xmsf02/xmw02  + 

5  lAr*(phAr+ahAr*log(t2) )*runiv*xmsfAr/xmwAr  + 

6  IH*  (phH+ahH*log(  t2  )  )  *runiv*xmsfH/xinwH  + 

7  lOH* ( phOH+ahOH*log ( t2 ) +bhOH*t2+chOH*t2*t2/2 . +dhOH* 

8  t2**3  ./3  • +ehOH*t2**4  ./4  .  )  *J^univ*xmsfOH/xmwOH  + 

9  lH2*(phH2+ahH2*log(t2)+bhH2*t2+chH2*t2*t2/2.+dhH2* 

9  t2**3  ./3  . +ehH2*t2**4  ./4  •  )  ■*i^univ*xmsfH2/xmwH2  + 

1  INO* ( phNO+ahNO*log ( t2 ) +bhNO*t2+chNO*t2*t2/2 . +dhNO* 

2  t2**3  ./3  • +ehNO*t2**4  ./4  •  )  *i^univ*xmsf  NO/xmwNO  + 

3  10* ( phO+ahO*log ( t2 ) +bhO*t2+chO*t2* t2/2 . +dhO* 

4  t2**3  ./3  • +ehO*t2**4  ./4  .  )  *i^’Jniv*xmsfO/xmwO  + 

5  1H20* ( phH20+ahH20*log ( t2 ) +bhH20*t2+chH20*t2*t2/2 . 

6  +dhH20*t2**3 ./3 .+ehH20*t2**4 ./4 • ) 

7  *runiv*xinsf  H20/xmwH20 
end  if 

f I=phimix2-phimix2g 

slope= ( phimix2g -phimixl ) / ( t2 - tl ) 

t3=t2+f 1/slope 

if (abs( (t3-t2)/t2) .It.rmtol)  go  to  20 

tl=t2 

t2“t3 

phimixl=phimix2g 
iter=iter+l 
100  continue 

v/rite ( 6  ,*)' tempera ture  for  region  after  shock  (derived 
2  from  entropy)  failed  to  converge.  Iters=',tter 
stop 

20  continue 
t=t3 
return 
end 

c=================end  of  subroutine  entropy================^ 
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subroutine  enthgam( t , hmix, gmix) 
c 

c  computes  enthalpy  and  gamma  as  a  function  of 

c  temperature  based  on  coefficients  from  CET85  program 

c  (JANAF  tables  curve  fitted) 
c 

common/coefficient/ 

1  alN2,blN2,clN2,dlN2,elN2,plN2,hlN2, 

2  ahN2 , bhN2 , chN2 , dhN2 , ehN2 , phN2 , hhN2 , 

3  al02 , bl02 , cl02 , dl02 , el02 , pl02 , hl02 , 

4  ah02 , bh02 , ch02 , dh02 , eh02 , ph02 , hh02 , 

5  alAr , plAr , hlAr , 

6  ahAr, phAr , hhAr , 

7  alH,blH,clH,dlH,elH,plH,hlH, 

8  ahH,bhH,chH,dhH,ehH,phH,hhH, 

9  alOH,blOH,clOH,dlOH,elOH,plOH,hlOH, 

9  ahOH , bhOH , chOH , dhOH , ehOH , phOH , hhOH , 

1  alH2,blH2,clH2,dlH2,elH2,plH2,hlH2, 

2  ahH2 , bhH2 , chH2 , dhH2 , ehH2 , phH2 , hhH2 , 

3  alNO,blNO,clNO,dlNO,elNO,plNO,hlNO, 

4  ahNO, bhNO, chNO,dhNO, ehNO, phNO , hhNO, 

5  aiO,blO,clO,dlO,elO,plO,hlO, 

6  ahO, bhO, chO, dhO, ehO, phO, hhO, 

7  alH20,blH20,clH20,dlH20,elH20,plH20,hlH20, 

8  ahH20,bhH20,chH20,dhH20,ehH20,phH20,hhH20, 

9  lN2,102,lAr,lH,10H,lH2,lN0,10,lH20 

common/species/ 

1  xmsfN2 , xmwN2 , xmsf02 , xmw02 , xmsf Ar , xmwAr , 

2  xmsf  H,  xmwH,  xmsf  OH,  xmv/OH,  xmsf  H2  ,  xmwH2 ,  xmsf  NO,  xmwN’O, 

3  xmsfO, xmwO, xmsf H20, xmwH20, 

4  runiv,rmix 

c . - . - . 

c  h  and  cp  coefficients  from  cet85 
c 

if  ( t .  le  .  1000  .  )  tiien 

hcoef N2=hlN2+alN2*t+blN2*t*t/2 . +clN2*t**3 . /3 . +dlN2* 

2  t**4 ./4 . +elN2*t**5 ./5 . 

hcoef02  =  hl02+al02*t-t-bl02*t*t/2  .  +cl02*t**3  .  /3  .  +dl02* 

2  t**4./4.+el02*t**5./5. 

hcoef  A_'=hlAr-5-alAr*t 
hcoef H=hlH+alH*t 

hcoefOH=hlOH+alOH*t+blOH*t*t/2  .  +clOH*t**3  .  /3  .  +dlOIi* 

2  t**4 ./4 . +elOH*t**5 ./5 . 

hcoefH2=hlH2+alH2*t+blH2’^t*t/2  .+clH2*t**3  ./3  .+dlH2* 

2  t**4 ./4 . +elH2*t**5 ./5 . 

hcoef NO=hlNO+alNO*t+blNO*t*t/2 . +clNO*t**3 . /3 . +dlNO* 

2  t**4./4.+elNO*t**5./5. 

hcoefO=hlO+alO*t+blO*t*t/2 . +clO*t**3 . /3 . +dlO* 
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t**4 ./4 . +eiO*t**5 ./5 ■ 

hcoef H20=hlH20+alH20*t+blH20*t*t/2 . +clH20*t**3 . /3 . 
+dlH20*t**4 ./4 .+elH20*t**5 ./5 • 


cpcoefN2=alN2+blN2*t+clN2*t*^t+dlN2*t**3  .  +elN2*t**4  . 
cpcoef02=al02+bl02*t+cl02*t^t+dl02*t**3 . +el02*t**4 . 
cpcoef Ar=alAr 
cpcoef H=alH 

cpcoefOH=alOH+blOH*t+clOH*t*t+dlOH*t**3 . +elOH*t**4  . 
cpcoef H2=alH2+blH2*t+clH2*t*t+dlH2*t**3 . +elH2*t**4  . 
cpcoefN0=alN0+blN0*t+clN0*t*t+dlN0*t**3 . +elN0*t**4 . 
cpcoef0=al0+bl0*t+cl0*t*t+dl0*t**3 .+el0*t**4 . 
cpcoef H20=alH20+blH20*t+clH20*t*t+dlH20*t**3 . 

2  +elH20*t**4. 

c 

else 

c 

hcoef  N2=hhN2-fahN2*t+bhN2*t*t/2  .  +chN2*t**3  .  /3  .  +dhN2* 
2  t**4 ./4 . +ehN2*t**5 ./5 . 

hcoef 02=hh02+ah02*t+bh02*t*t/2 . +ch02*t**3 . /3 . +dh02* 
2  t**4 ./4 ■ +eh02*t**5 ./5 . 

hcoef Ar=hhAr+ahAr*t 
hcoef H=hhH+ahH*t 

hcoef 0H=hh0H+ah0H*t+bh0H*t*t/2 . +chOH*t**3 ./3 . +dhOH* 
2  t**4./4.+ehOH*t**5./5. 

hcoef H2=hhH2+ahH2*t+bhH2*t*t/2 . +chH2*t**3 . /3 . +dhH2* 
2  t**4./4.+ehH2*t**5./5. 

hcoef N0“hhN0+ahN0*t+bhN0*t*t/2 . +chN0*t**3 ./3 . +dhNO* 
2  t**4./4.+ehNO*t**5./5. 

hcoef 0=hh0+ah0*t+bh0*t*t/2 . +ch0*t**3 . /3 . +dhO* 

2  t**4  . /4  -  i-eh0*t**5  ./5  - 

hcoef  H20=hhH20+ahH20*t+bhH20*t*t/2  .■i-chH20*t**3  ./3  . 

2  +dhH20*t**4 ./4 . +ehH20*t**5 . /5 . 

c 

cpcoefN2=ahN2+bhN2*t+chN2*t*t+dhN2*t**3 .  +ehN’2*t**4  . 
cpcoef02=ah02+bh02*tfch02*t*t  {■dh02^t**3  .  +eh02*t**4  . 
cpcoef Ar=ahAr 
cpcoef H=ahH 

cpcoefOH=ahOH+bhOH*t+chOH*t*t^dhOH*t**3  .  +ehOH'^t**4  . 
cpcoef H2=ahH2+bhH2*t+chH2*t*t+dhH2*t**3 . +ehH2*t**4  . 
cpcoefN0=ahN0+bhN0*t+chN0*t*t+dhN0*t**3 . +ehN0»t**4 . 
cpcoef0=ah0+bh0*t+ch0*t*t+dh0*t**3 .+eh0*t**4 . 
cpcoefH20=ahH20+bhH20*t  +  chH20'»t*t+dhH20*t**3  . 

2  +ehH20*t**4. 

c 

end  if 
c 

c  Compute  h  for  each  species.  Molar  values  of  h  (t=0K) 
c  must  be  found  for  each  species.  These  can  be  found 


2 

2 

c 

c 
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c 

c 


in  the  JANAF  tables . 


c 

c 


hN2=hcoef N2*runiv+8669000 . 
h02=hcoef02*runiv+8682000 . 
hAr=hcoef Ar*runiv+6196450 . 
hH=hcoef H*runiv+6196504 . 
hOH=hco€fOH*runiv+8815688 . 
hH2=hcoef H2*runiv+8468416 . 
hNO=hcoef NO*runiv+9192248 . 
hO=hcoefO*runiv+6727872 . 
hH20=hcoef H20*runiv+9904000 . 

Compute  h  of  mixture 

hmix=lN2*hN2*xmsf  N2/xmv/N2-i- 

2  102*h02*xmsf02/xmv/02  + 

3  lAr*hAr*xmsf Ar/xmwAr+ 

4  lH*hH*xmsf H/xmwH+ 

5  10H*h0H*xmsf0H/xmv/0H+ 

6  lH2*hH2*xmsfH2/xmwH2-r 

7  lNO*hNO*xmsf NO/xmwNO+ 

8  10*h0*xmsf0/xmw0-5- 

9  lH20*hH20*xmsfH20/xmwH20 

Compute  cp  of  mixture 

opmix=lN2*cpcoefN2*runiv*xmsfN2/xmv/N2 

2  +I02*cpcoef02*runiv*xmsf02/xmw02 

3  +lAr*cpcoef Ar*runiv*xmsf Ar/xmwAr 

4  +lH*cpcoef H*runiv*xmsf H/xmwH 

5  +iOH*cpcoefOH*runa  v*xmsfOH/xmv/OH 

6  +lH2*cpcoef  H2*runiv*xmsf  H2/xmv/H2 

7  +lNO*cpcoef  NO*runiv*xmsfNO/xmv/NO 

8  +10*cpcoef0*runiv*xmsf0/xmv/0 

9  +lH20*cpcoef H20*runiv*xmsf H20/xmwH20 

Compute  cv  and  gamma  of  the  mixture 

cvmix=cpmix- rmix 
gmix=cpjnix/cvmix 

return 

end 

===============  end  of  subroutine  enthgam==== 
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4^u)to 


subroutine  imperfect  (j,g) 
c 

c  computes  the  primitive  variables  from  the  e 
c  vector  based  on  a  thermally  perfect  gas 
c 

parameter  (jdim=501) 
common/species/ 

xmsf  N2 ,  xmwN2 ,  xms  f 02  ,  xmw02 ,  xms  f  Ar ,  xmv/Ar , 

2  xms f H, xmwH, xmsf OH, xmwOH, xmsf H2, xmwH2 , xmsf NO, xmwNO, 

3  xms f 0,  xmwO,  xmsf  H20,  xmv/H20, 

4  runiv,rmix 

common/solution/  cn , olddx, pamb, betafd, rmtol , xstop, 

phi , x( jdim, 2 ) , y ( jdim, 2 ) , rh( jdim, 2),u(jdim,2), 
v( jdim, 2 ) , p( jdim, 2 ) , t ( jdim, 2 ) , rhe ( jdim, 2 ) , 
xmach( jdim, 2 ) 

common/riem/  rhO ( jdim) , rh2 ( jdim) , rh4 (jdim) , rh6 (jdim) , 
uO( jdim) ,u2( jdim) ,u4( jdim) ,u6( jdim) , xmach2 ( jdim) , 
vO( jdim) ,v2( jdim) ,v4( jdim) , v6( jdim) , xmach4 ( jdim) , 
pO( jdim) ,p2( jdim) ,p4( jdim) ,p6( jdim) , 
rheO( jdim) ,rhe2( jdim) ,rhe4( jdim) ,rhe6( jdim) , 
slplO( jdim) ,slp22( jdim) ,slp36( jdim) 

common/e  f  vector/ 

1  el( jdim,2) ,e2(jdim,2),e3(jdim,2) ,e4( jdim,2) , 

2  f 1 ( jdim, 2 ) , f 2 ( jdim, 2 ) , f 3 ( jdim, 2 ) , f 4 ( jdim, 2) 

common/compf Igl/ 

1  jmarch, ibnd, jattach, ioptim, itermax, iteropt, 

2  imax,  jnozz ,  ix ,  knozzle ,  nv/pnts ,  kucowl ,  klcowl , 

3  ipsarc, igrid, 

4  islpnoz, jucowi, jlcowl, ipackin, ipackex, 

5  jpackin , jpackex , 

6  iriem, irmbnd, inpnoz, ivint, ivext, nixint, 

7  nixext,kbot, 

5  method, limit, ipastuc, ipastlc, ixcowl, irmfix, 

9  itri.em,  itmoc, 

9  inpcwl , icompar , ktrk , idoext 
c 

c . - . 

c  j  =  l 

iter=0 

itmax=50 

if  (  j  .eq  .  jnozz  .  or  .  j  .eq  .  j  icov/1 )  then 

tl=pO( j)/rhO( j)/rmix 

else 

tl=p6( j )/rh6 ( j )/rmix 
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endif 

t2=tl+5a. 

call  enthgam  (tl/hl/gl) 

ua= (e2 ( j , 2)+sqrt (e2 ( j , 2) *e2 ( j , 2 ) -4 . *el ( j / 2 ) * 

2  el  ( j  ,  2 )  *rinix*tl)  )/2  ./el{  j  ,  2 ) 
fa=hl*el (j,2)+.5*el(j,2) *ua*ua+ . 5*e3 ( j / 2 ) * 

2  e3(j,2)/el( j,2) 

10  call  enthgam  (t2,h2/g2) 

ub= (e2 ( j , 2)+sqrt (e2 ( j , 2) *e2 ( j , 2 ) -4 . *el ( j / 2 ) * 

2  el( j , 2 ) *rmix*t2 ) )/2 ./el( j / 2 ) 
fb=h2*el( j , 2 )+ . 5*el( j , 2)*ub*ub+ . 5*e3 ( j , 2) * 

2  e3( j,2)/el( j,2) 

slope=(fb-fa)/(t2-tl) 
t3=t2+(e4 ( j , 2) -fb) /slope 
if (abs( (t3-t2)/t2) .le.rmtol)  go  to  20 
if ( iter . gt . itmax)  then 

write ( 6 , * ) ' sub .  imperfect  failed  to  converge,  plane 
2  ix= ' , ix 

write ( 6 , * ) ' and  stops  at  node  j=',j 
stop 
endif 
fa=fb 
tl=t2 
t2=t3 

iter=iter+l 
go  to  10 
20  continue 
g=g2 

U(j,2)=(e2(j,2) +sqrt {e2{j,2)*e2(j,2)-4.*el(j,2)* 

2  el ( j , 2 ) *rmix*t3 ) )/2 . /el ( j , 2 ) 

v( j,2)=e3( j,2)/el( j,2) 
if (abs (v(j,2)).lt.l.e-3)  v(j,2)=0.0 
if (abs (e3 (j,2)).lt.l.e-3)  e3(j,2)=0.0 
p(j,2)=e2( j,2)-el(j,2)*u( j,2) 
rh(j,2)=el(j,2)/u( j,2) 

rhe (j,2)=rh{j,2) *h2 -p ( j , 2 ) + ( rh ( j , 2 ) /2 . ) * 

2  (u(j,2)*u(j,2)+v( j,2)*v(j,2)) 

ax=asnd(g2,p( j , 2) ,rh( j,2) ) 
vmag=sqrt (u(j,2)*u{j,2)+v(j,2)*v(j,2)) 
xmach( j , 2)=vmag/ax 
t(j,2)=t3 

if ( ix . eq . 39 . and. ( j . eq . 50 .or . j . eq . 51 ) ) then 
endif 

return 

end 
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subroutine  enthgam(t,hmix,gmix) 
c 

c  computes  enthalpy  and  gamma  as  a  function  of 

c  temperature  based  on  coefficients  from  CET85  program 

c  (JANAF  tables  curve  fitted) 

c 

common/coefficient/ 

1  alN2 , blN2 , clN2 , dlN2 , elN2 , plN2 , hlN2 , 

2  ahN2 , bhN2 , chN2 , dhN2 , ehN2 , phN2 , hhN2 , 

3  al02 , bl02 , cl02 , dl02 , el02 , pl02 , hl02 , 

4  ah02 , bh02 , ch02 , dh02 , eh02 , ph02 , hh02 , 

5  alAr, plAr, hlAr , 

6  ahAr , phAr , hhAr , 

7  lN2,102,lAr 
common/species/ 

xms  f N2 , xmwN2 , xms  f 02 , xmw02 , xmsf Ar , xmwAr , 

2  runiv,rmix 

c . . - . 

G  h  and  cp  coefficients  from  cet85 
c 

if ( t . le . 1000 . )  then 

hcoef N2=hlN2+alN2*t+blN2*t*t/2 . +clN2*t**3 . /3 . +dlN2* 

2  t**4./4.+elN2*t**5./5. 

hcoef 02=hl02+al02*t+bl02*t*t/2 . +cl02*t**3 . /3 . +dl02* 

2  t**4./4.+el02*t**5./5. 

hcoefAr=hlAr+alAr*t 

cpcoef N2=alN2+blN2*t+clN2*t*t+dlN2*t**3 . +elN2*t**4 . 
cpcoef02=al02+bl02*t+cl02*t*t+dl02*t**3 . +el02*t**4 . 
cpcoef Ar=alAr 
else 

hcoef N2=hhN2+ahN2*t+bhN2*t*t/2 . +chN2*t**3 ./3 . +dhN2* 

2  t**4 . /4 . +ehN2*t**5 ./5 . 

hcoef02=hh02+ah02*t+bh02*t*t/2 . +ch02*t**3 ./3 . +dh02* 

2  t**4 ./4 .+eh02*t**5 ./5 . 

hcoef Ar=hhAr+ahAr*t 

cpcoef N2=ahN2+bhN2*t+chN2*t*t+dhN2*t**3  .  +ehN2’*'t**4  . 
cpcoef 02=ah02+bh02*t+ch02*t*t+dh02*t**3 . +eh02*t**4 . 
cpcoef Ar=ahAr 
end  if 
c 

c  Compute  h  for  each  species.  Molar  values  of  h  (t=298K) 
c  must  be  found  for  each  species.  Naturally  occuring 
c  species  will  not  have  heats  of  formation, 
c 

hN2=hcoef N2*runiv+8669000 . 
h02=hcoef02*runiv+8682000 . 
hAr=hcoef Ar*runiv+6196000 .45 
c 

c  Compute  h  of  mixture 
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c 

hinix=lN2*hN2*xinsfN2/xmwN2+ 

2  I02*h02*xmsf02/xmw02+ 

3  lAr*hAr*xmsf Ar/xmwAr 
c 

c  Compute  cp  of  mixture 
c 

cpmix=lN2*cpcoefN2*runiv*xmsfN2/xmwN2 

2  +102*cpcoef02*runiv*xmsf02/xmw02 

3  +lAr*cpcoef Ar*runiv*xmsf Ar/xmwAr 
c 

c  Compute  cv  and  gamma  of  the  mixture 
c 

cvmix=cpmix - rmix 
gmix=cpmix/cvmix 

return 

end 

c==================-===============--===================== 

subroutine  tfromh(h2, tt, t2) 
c 

c  back  calculate  temperature  given  enthalpy  using  the 
c  Newton -Raphson  iteration  method 
c 

common/coefficient/ 

1  alN2,blN2,clN2,dlN2,eli':2,plN2,hlN2, 

2  ahN2 , bhN2 , chN2 , dhN2 , ehN2 , phM2 , hhN2 , 

3  al02 , bl02 , cl02 , dl02 , el02 , pl02 , hl02 , 

4  ah02 , bh02 , ch02 , dh02 , eh02 , ph02 , hh02 , 

5  alAr, plAr , hlAr, 

6  ahAr , phAr , hhAr , 

7  lN2,102,lAr 

common/species/ 

xmsf N2 , xmwN2 , xmsf02 , xmw02, xmsf Ar , xmwAr , 

2  runiv,rmix 

c . - . 

itmax=300 

tol=l.e-6 

do  100  iter=l,itmax 
if ( tt . le . 1000 . )  then 

hcoefN2=hlN2+alN2*tt+blN2*tt*tt/2 .+clN2*tt**3 ./3 • 
2  +dlN2*tt**4./4 .+elN2*tt**5./5. 

hcoef02=hl02+al02*tt+bl02*tt*tt/2 . +cl02*tt**3 . /3 . 
2  +dl02*tt**4 ./4 . +el02*tt**5 ./5 . 

hcoef Ar=hlAr+alAr*tt 
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else 

hco;fN2=hhN2+ahN2*tt+bhN2*tt*tt/2 . +chN2*tt**3 ./3 . 

2  +d[hN2*tt**4./4.-*-ehN2*tt**5./5. 

hcoef02=hh02+ah02*tt+bh02*tt*tt/2 . +ch02*tt**3 . /3 - 
2  +dh02*tt**4./4.+eh02*tt**5./5. 

hcoefAr=hhAr+ahAr*tt 
end  if 

if ( tt . le . 1000 . )  then 

dhcoef  N2=alN2+blN2*tt+clN2*tt**2  .  /2  .  +c:iN2* 

2  tt**3 ./3 .+elN2*tt**4 ./4 . 

dhcoef02=al02+bl02*tt+cl02*tt**2 . /2 . +dl02* 

2  tt**3 ./3 .+el02*tt**4 ./4 . 

dhcoefAr=alAr 
else 

dhcoefN2=ahN2+bhN2*tt-i-chN2*tt**2  ./2  •  +dhN2* 

2  tt**3 ./3 .+ehN2*ct**4 ./4 . 

dhcoef02=ah02+bh02*tt+ch02*tt**2 ./2 . +dh02* 

2  tt**3 ./3 .+ehC2*tt**4 ./4 ■ 

dheoef Ar=ahAr 
end  if 

hN2=hcoefN2*runiv+8669000 . 
h02=hcoef02*runiv+8682000 . 
hAr=hcoefAr*runiv+6196450 . 

hinix=lN2*hN2*xmsfN2/xmwN2+ 

2  102*h02*xmsf02/xmw02+ 

3  lAr*hAr*xmsfAr/xmwAr 

ddhN2=dhcoefN2*runiv 
ddh02=dhcoef02*runiv 
ddhAr=dhcoef Ar*runiv 

dhmix=-  ( lN2*ddhN2*xmsfN2/xmv/N2+ 

2  102*ddh02*xmsf02/xmw02+ 

3  lAr*ddhAr*xmsfAr/xmv;Ar ) 

fh=h2-hmix 

t2=tt-fh/dhmix 

c  write(6,*) ' t2  iterated  from  tfrmoh',t2 

if  (abs( (tt-t2)/t2) . It. tol)  go  to  20 
tt=t2 

iter=iter+l 
hmix=0 . 

100  continue 

write( 6,*) 'temperature  from  enthalpy  failed  to 
2  converge' 

stop 
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*###  n  it  t  ###  fi  ##i  sp  if  #  n##  s  if  ######## 

*  WARNING:  THIS  PROGRAM  IS  VERY  SENSITIVE  TO  INITIAL 

*  GUESSES  OF  EPSILON.  SHOULD  BE  SET  CLOSE  TO 

*  ACTUAL  VALUE ! 

*  if  it  ft  it  #  if  if  if  if  it  it  #  it  ##  if  fi  #  if  if  if  if  #  If  if  if  ft  t;  it  #  if  #  ft  fi  if  ft  ft  if  it  i  it  if  if  it  if  it  ft  it  if  if  fi  if  #  if  if  if  ii  ft 

*****-********it********************^******it*************'**** 

*  Thi.v  program  will  solve  the  "exact"  oblique  shock  * 

*  problem  for  an  imperfect  gas  (air),  patterned  after  the* 

*  example  problem  in  Zucrow  and  Hoffman  example  7-10.  * 

*  The  output  of  this  is  used  to  form  the  IV  line  that  * 

*  is  used  to  validate  the  imperfect  gas  model.  * 

************<r<i***#******'*<:*******l<:***1lc****'*************1t*:fc* 

program  main 
common/coef  f ic ient/ 

1  a 1N2 , blN2 , clN2 , dlN2 , elN2 , plN2 , hlN2 , 

2  ahN2 , bhN2 , chN2 , dhN2 , ehM2 , phN2 , hhN2 , 

3  ai02 , bl02 , cl02 , dl02 , el02 , pl02 , hl02 , 

4  ah02 , bh02 , ch02 , dh02 , eh02 , ph02 , hh02 , 

5  alAr , plAr, hlAr , 

6  ahAr,phAr,hhAr, 

7  lN2,102,lAr 

common/species/xmsf  N2 ,  xn;wN2 ,  xmsf  02 ,  xmv/02 ,  xmsf  Ar ,  xmwAr , 
2  runiv,rmix 

c . 

tol-l.e-10 

itmax=200 

1N2=1 

102=1 

lAr=l 

c 

c  xmw  -  molecular  weight 
c  xmsf-  mass  fraction 

c  a,b,c  etc.  CET  coef f icicents  (1-  low  temp,  h-  high  temp) 
xmwN2=28.013 
xmsfN2=. 75556737622 
alN2=3. 7044177 
blN2=- .142187536-2 
clN2=.2.8670392e-5 
dlN2=- . 12028b85e-8 
elN2=- .139540776-13 
plN2=2. 2336285 
hln2=- . 10640795e4 

ahN2=2. 8532899 
bhN2= . 16022128e-2 
chN2=- . 62936893e-6 
dhN2=.11441022e-9 
ehN2=- .780574656-14 
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phN2=6. 3964897 
hhN’2=-  .  89008093e3 

xmw02=32 . 
xinsf02=.  231605141 
al02=3 .7837135 
bl02=- .302336346-2 
cl02- . 994927516-5 
dl02=- .981891016-8 
6102=. 330318256-11 
pl02=3. 6416345 
hl02=- .1063810764 

ah02=3 .6122139 
bh02=. 748531666-3 
ch02=- .198206476-6 
dh02= . 337490086-10 
6h02=- .239073746-14 
ph02=3 .6703307 
hh02=- . 1197815164 

xmwAr=39 . 944 
xinsfAr=. 01282748278 
alAr=2.5 
plAr=4 . 3660006 
hlAr=- .7453749863 

ahAr=2 . 5 
phAr=4 . 3660006 
hhAr=- .7453750263 

rmix=287 . 09 
runiv=8314 . 34 
c 
c 

pi=3 . 14159265359 
degrad=pi/180 . 

*  initial  flow  properti6s 

write(6,*) 'input  delta  in  degrees' 
read{5,*)  delta 
delta=deita*d6grad 

write ( 6 ,*)' input  temperature  in  Kelvin' 
read (5,*)  tl 

write ( 6, *)' input  mach  number' 
read (5,*)  xml 

c  write(6,*) 'input  velocity  (m/s)' 

c  read(5,*)  vl 

write(6,*} '  input  pressure  (N/m''2)' 
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read(5/*)  pi 

c  .  write(6,  *) 'input  density 

c  read(5,*)  rhol 

wr.Lte(  6 ,  * ) '  input  first  initail  guess  for  epsilon' 

read{5..>f)  epsl 

epsl=epsl*degrad 

write( 6 ,  ' input  second  initail  guess  for  epsilon' 

read(5,-0  eps2 

eps2=eps2*deg.rad 

★ 

rhol=pl/tl/rmix 

call  enthgam( tl/ hmix, gmix) 

al=sqrt ( gmix*rmix* t 1 ) 

vl=xinl*al 

hl=hmix 

write ( 6, *) 'delta  =  ', delta/degrad 

write(  6 ,  *) 'Ml  =  ',xinl 

write( 6, * ) ' vlmag  =  ',vl 

write(6, *) 'al  =  '.al 

write(6,  *) 'pi  ',pl 

write(6/ *) ' tl  =  ',tl 

write(6, *) 'rhol  =  ',rhol 

write ( 6 ,  * ) ' gaimRal= ' ,  gmix 

write(6,*)' ' 

do  10  iter=l,itmax 

fepsl=(  (gmix+1  .)/2  .  *xnil’*xml/(xml*xml  -  sin (epsl) 

2  *sin(epsl) -1. ) ■!• )*tan(epsl) 

feps2=( (gmix+1. )/2 . *xml*xml/(xml*xml*sin(eps2) 
*.9in(eps2)  -1 . )  - 1 . )  *tan  (  eps2 ) 
slope=  ( fep,s2  -  fepsl)/(eps2-epsl) 
eps3=eps2+( l./tan (delta) - feps2)/slope 
c  write( 6 ,») 'epsilon  iterations ', eps3 

if  (abs ( (eps3 -eps2 )/eps2 ) . It . tol )  go  to  20 
epsl=eps2 
eps2=eps3 
■ ter=iter+l 

i  (iter.gt.itmax)  then 

V..  ■  te(6,  *) 'epsilon  did  not  converge' 

st 

encif 

10  continue 
20  epsl=eps3 
c 

50  xnml=xml*sin(epsl ) 

xnm2=sqrt ( (yi:ml*xnml+2  ./(gmix  - 1 . ) )/( 2  .  *gmix/(gmix- 1 .  ) 
2  *xnml*xnml-l . ) ) 
xm2g=xnm2/(sin(epsl'deita) ) 

rhorat=( gmix+1. ) *xnml*xnml/( 2 . + (gmix- 1 . )*xnml*xnml) 
rho21=rhol*rhorat 
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trat=(2.*gmix/(gmix+l. )  *xnnil*xnml-  ( (gmix-1.  )/gmix-M.  ) ) 
2  *( ( (gmix-1 . )/(gmix+l . ) )+2 ./( (gmix+1 . ) *xnml*xnml) ) 

t2g=trat*tl 

prat=2*gmix/(gmiXTl . )  *xnml*xnml-  { (ginix-1 .  )/(gmix+l .  ) ) 
p2g=prat*pl 

vnl=vl*sin (epsl) 
vt=vl*cos (epsl ) 

p2=pl+rhol*vnl*vnl* (1 . -rhol/rho21 ) 
h2=hl+vnl*vnl/2 .*(!.-( rhol/rho21 ) * ( rhol/rho21 ) ) 
call  tfromh (h2 , t2g, t2 ) 
rho22=p2/rmix/t2 

if  (abs( (rhc22-rho21)/rho21) . It . tol)  go  to  40 
p2=pl+rhol*vnl*vnl* ( 1 . -rhol/rho22 ) 
h2=hl+vnl*vnl/2 . ^ ( 1 . - ( rhol/rho22 ) * ( rhol/rho22 ) ) 
call  tfromh(h2 , t2g, t2) 

c  Secant  method  for  iteratively  finding  density 
30  rho23=p2/rmix/t2 

if  (abs( (rho23-rho22)/rho22) . It. tol)  go  to  40 

slope=(rho22-rho23 )/(rho21-rho22) 

rho2c=rho23+slope*(rho2c.  -  rho22) 

if  (abs(  (rho2c-rho23)/rho.23)  .It.tol)  go  to  40 

rho21=rho22 

r«io22=rho23 

p2=pl+rhol*vnl*vnl* ( 1 . -rhol/rho2c) 
c  write (6, *)' iterations  of  p'/p2 

h2=hl+vnl*vnl/2 . * ( 1 . - ( rhol/rho2c ) * ( rhol/rho2c ) ) 
call  tfromh(h2, t2g, t2) 
go  to  30 

40  continue 
rho2=rho2c 
vn2=rhol/rho2*vnl 
v2=sqrt ( vt*vt+vn2*vn2 ) 
call  enthgam(t2,hmix,gmix) 
a2=sqrt ( gmix*rmix*t2 ) 
xm2=v2/a2 

c  write(6,*) 'mach2=' ,xm2 
eps2=delta+asin ( vn2/v2 ) 

if  (abs( (eps2-epsl)/epsl) .It. tol)  go  to  60 
epsl=eps2 
go  to  50 
60  continue 

write (6,*) 'epsilon  =' ,eps2/degrad 
write{6, *) 'M2  =',xm2 
write(6, *) ' v2mag  =',v2 
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